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RELATIONSHIP  BETWEEN  NEAR-FIELD  AND  TELESEISMIC 
OBSERVATIONS  OF  SEISMIC  SOURCE  PARAMETERS 


Brief  Summary 


I.  Objective 

The  overall  objective  of  this  investigation  has  been  to  establish  the 
connection  between  source  mechanism  parameters  inferred  from  seismic 
measurements  made  relatively  near  to  earthquakes  and  explosions  and  those 
inferred  from  teleseismic  observations. 

II . Major  Accomplishments 

The  following  brief  summary  outlines  the  accomplishments  under  this 
grant.  We  have: 

1.  found  that  there  are  significant  variations  in  body-wave  and  surface 
wave  excitation  for  explosions  of  comparable  magnitude  detonated  in  similar 
source  media  and  located  only  a few  kilometers  apart  at  NTS  (Yucca  Flats  and 
Pahute  Mesa).  These  observed  differences  in  excitation  and  energy  part- 
itioning imply  either  a rapidly  varying  tectonic  stress  field  in  the  source 
region  and/or  strongly  varying  patterns  of  near-source  fracturing.  Source- 
generated P-wave  coda  of  two  minutes  duration  or  more  are  observed  at  NORSAR 
to  vary  significantly  in  strength  and  character  among  nearby  Yucca  and  Pahute 
events,  possibly  due  to  short-period  surface  wave  to  P-wave  scattering; 

2.  shown  that  relatively  few  SRO-type  stations  can  be  used  to  compare 
source  mechanisms  for  suites  of  events  of  varying  size  in  a given  source  region. 


For  example,  foreshocks  and  aftershocks  of  the  Utah-Idaho  border  and  Oroville, 
California  sequences  of  1975,  as  recorded  at  the  Albuquerque  SRO  station,  have 
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been  analyzed  and  compared;  In  both  cases  the  principal  foreshock  exhibited 
the  same  mechanism  as  the  main  shock,  while  the  aftershocks  are  more  varied 
in  mechanism.  For  both  sequences  we  were  able  to  match  the  observed  Love  and 
Rayleigh  waveforms  and  spectra  as  the  mechanism  changed.  Average  path 
dispersion  (hence  structure)  is  a useful  by-product  of  the  analysis; 

3.  developed  a promising  approach  for  using  the  entire  short-period 
signature  for  identifying  short-time- lag  double  events;  the  technique  will 
work  best  at  regional  distances  where  relatively  broad-band  surface  wave 
(e.g.  Lg,  Rg)  as  well  as  body  waves  (Pn,  Pg,  Sg)  signals  are  recorded; 

4.  developed  further  the  surface-wave  inversion  techniques  we  devised 
for  inferring  lateral  heterogeneity  in  structure  using  any  combination  if  Love 
and  Rayleigh  group  velocity  or  phase  velocity  dispersion  for  as  many  modes 

as  are  observed; 

5.  developed  a finite  element  code  that  can  compute  theoretical  wave- 
forms that  would  be  observed  at  any  point  of  interest  from  an  arbitrary 
source  extended  in  space  and  time  in  geologically  heterogeneous  structures. 

In  addition  we  have  been  able  to  combine  this  code  with  the  Haskell  method 
to  propagate  waves  from  a simple,  plane-layered  region  into  and  across  a 
complicated  structure.  Using  this  code  the  1971  San  Fernando  earthquake 
was  simulated  using  static  and  dynamic  fault  models  with  a variable  rupture 
velocity.  We  obtained  good  agreement  with  observed  near-field  observations 
for  a rupture  velocity  of  about  2 km/sec.  We  have  obtained  digitized  records 
for  additional  near-field  sites  and  are  now  attempting  to  interpret  them 
with  the  aid  of  the  finite  element  calculations;  this  involves  interpolation 
of  the  digital  data  to  a uniform  sampling  rate  in  time  and  filtering  to 
produce  a wave-form  with  the  same  spectral  band  as  the  numerical  code 
calculations; 
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6.  far-field  surface  wave  amplitude  spectra  have  been  studied  for  their 
use  in  determining  seismic  source  parameters.  Using  theoretical  spectra 
generated  for  a double  couple  source  in  a plane  layered  half-space,  the 
minimum  residual  fit  to  the  observed  spectra  was  obtained,  yielding  source 
parameter  estimates;  and 

7.  developed  a software  package  for  analysis  of  SRO  data  tapes.  In 
addition  to  standard  types  of  filtering,  matched  filtering,  spectral 
decomposition,  and  dispersion  analysis,  we  have  implemented  routines  for 
polarization  filtering  to  isolate  individual  phases  and  stacking  of  events. 

This  grant  has  supported  wholly  or  in  part  two  completed  Ph.D.  theses 
and  three  that  are  near  completion.  In  addition,  the  soft-ware  developed 
under  this  grant,  particularly  the  inversion  programs,  has  been  desemminated 
rather  widely  and  it  is  being  used  by  several  organizations  on  ARPA/AFOSR 
projects . 

III.  Pub  1 lent  ions 

Publications  based  on  the  work  summarized  above  are  listed  on  the 
following  two  pages. 
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Contained  in  this  Appendix  are  examples  to  Illustrate  some  of  the  results 
discussed  In  the  Summary  of  Accomplishments  but  not  Incorporated  Into  the 
references,  publications,  or  reports  that  appear  elsewhere  In  this  document. 


TABLE  A- I 

This  Table  summarizes  Che  variations  in  Love  and  Rayleigh  wave  magnitudes 
expected  as  a function  of  source  structure,  degth,  and  fault  type  (tectonic 
stress);  VSS  ■ vertical  strike-slip;  45DS  ■ 45  dip-slip;  periods  as  given 
and  amplitudes  were  measured  from  synthetic  seismograms  for  3000  km  distance. 
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TABU!  II 

BASIC  SRO  PROGRAMS : 


I  HEADFRDUMP. 

READS  MULTIPLE  FILE  SRO  TAPES  AND  PRINTS  HEADER  LABEL  FOR  EVERY 
RECORD. 

7 SPSUBSET. 

READS  SEOHENT  OF  AN  SRO  TAPE.  EXTRACTS  SPECIFIED  SHORT  PERIOD  WINDOW, 
DEMAONIFIES  DATA,  AND  PRODUCES  A STANDARD  FORMAT  SUBSET  TAPE  AND  A 
PLOT  TAPE. 

3 LPSUPSET. 

SIMILAR  TO  SPSURSET  WITH  ADDITIONAL  FEATURE  OF  DENULT I PI. FX I HQ  THREE 
COMPONENT  LONO  PERIOD  DATA  AND  ROTATING  HORIZONTAL  COMPONENTS  TO  0 
SPECIFIED  DIRECTION. 


SRO  DATA  PROCESSING  PROORAMS: 


1 COLLAPSE. 

AN  OLD  SDL  PROORAM  TO  NATCHFII.TER  LONG  PERIOD  DATA  USING  REAL 
OR  SYNTHETIC  SIONALS  AS  FILTFRS.  ADDITIONAL  FEATIIRFS  AI.I.OU 
COMPUTATION  OF  SMALL ‘EVENT  SPECTRA,  BEAM I NO  OK  ARRAY  DATA,  ETC. 

3 NARROUBANO • 

PERFORMS  NARROW  BAND  FILTERING  IN  THE  FREQUENCY  DOMAIN  FOR 
EXTRACTION  OF  OROUP  ARRIVAL  TIMES.  JUDICOIIS  CHOICE  nF  WIPFR 
PASS  BAND  ALLOWS  SEPERATION  OF  NORMALLY  AND  INVERSELY  DISPERSED 
ARRIVALS. 

3 SEISMOPRINT. 

PRODUCES  3-D  CONTOUR  PLOT  OF  POWER  VS  FREQUENCY  AS  A FUNCTION  OF 
TIME.  OPTIONALLY,  PROORAM  ALSO  PRODUCES  INSTRUMENT  CORRECTED  PLOTS 
OF  POWER  VS  FREQUENCY  AVERAGED  OVER  A GIVEN  TIME  WINDOW. 

4 LOO-LOOPLOT,  . 

COMPUTES  POWER  SPECTRUM  OF  A LENGTH  OF  SEISMIC  DATA  AND  PRODUCES  A 
CALCOMP  PLOT  OF  LOO  POWER  VS  LOO  FREQUENCY,  THUS  FACILITATING  THE 
DETERMINATION  OF  CORNER  FREQUF.NCY . 

3 LTOR. 

COMPUTES  THE  POWER  SPECTRAL  DENSITY  OF  THF  LOVE  AND  RAYLEH.lt  UA'.'F 
(VERTICAL  COMPONENT)  PORTION  OF  A LONG  PERIOD  SEISMOGRAM  AND  PRODUCES 
A PLOT  Of  THE  RATIO  L/R  AS  A FUNCTION  OF  PERIOD. 


ANCILLARY  PROORAMS. 


1 XWIP2M. 

INVERTS  OBSERVED  SURFACE  WAVE  DISPERSION  DATA  TO  DETERMINE  THE  BEST 
SINGLE  OR  COMPOSITE.  FLAT  OR  SPHERICAL,  EARTH  NCSEL.  INPUT  CONSISTS 
OP  A GUESS  MODEL,  SINGLE  OR  COMPOSITE.  THE  OBSERVED  DATA  AND  FS1JMATFS 
OF  THEIR  VARIANCES.  AND  A UEIOHTINO  MATRIX  EXPRESSING  CONFIDENCE  IN 
THE  TRIAL  MODEL.  DATA  CAN  BE  LOVE  OR  RAYLF.IGM  WAVE  PHASE  AND  GROUP 
VELOCITIES,  FUNDAMENTAL  OR  HIGHER  MODE. 

2 HARPING. 

COMPUTES  LOVE  AND  RAVLEIOH  WAVE  PHASE  AND  GROUP  VELOCITIES  AMO  THEIR 
PART7ALS  WITH  RESPECT  TO  THE  MODEL  PARAMETER*  FOR  A PLANf-PARAtl FI 
LAYERED  HALFSPACC. 

3 XTRTCH. 

CONVERTS  INPUT  FLAT  FARTN  MODEL  INTO  A CORRESPONDING  PSFUDn-SPMFPTCM 
MODEL  FOR  USE  WITH  HARPING. 

4 RAYPtINCH  AMft  L nUPPlINCM 

COMPUTES  THE  RAVI CIOH  OR  LOVE  WAVE  DISPLACEMENT- STRESS  OIIANTITJFS  AT 
AN  IMAGINARY  SURFACE.  CONTAININO  THE  POINT  SOURCE  FOR  A Pt.ANE-NARAl  l.EI. 
LAYERED  HALFXPACF.. 

3 SURXIM. 

COMPUTES  SYNTHETIC  LOVE  AMD  RAYLEIGH  WAVE  SEISMOGRAMS  FROM  EARTHQUAKE 
AND  EXPLOSIVE  SOURCES  IN  A PLANE-PARALLEL  LAYERED  EARTH  MODFI. 
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Figure  A-l.  Short-period  and  long-period  seismograms  recorded  at  the 

Albuquerque  SRO  station  for  two  nuclear  explosions  (Keelson 
and  Esrom)  located  within  approximately  2 km  of  one  another 
at  Yucca  Flats  and  detonated  about  20  minutes  apart  in  time. 
Note  that  the  levels  of  excitation  for  the  short-period 
phases  (Pn,Pg,Lg)  are  approximately  the  same  for  the  two 
events,  whereas  there  is  a very  large  difference  in  both 
Love  and  Rayleigh  wave  excitation  between  the  two  with  the 
second  (Esrom)  exhibiting  the  large  surface  waves  (tectonic 
release).  Similar  results  were  obtained  for  other  (SDCS) 
stations . 


KEELSON 


Figure  A-l-c 
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sponsorship  from  this  grant. 


ABSTRACT 


DYNAMIC  FINITE  ELEMENT  ANALYSIS  WITH  APPLICATIONS 
TO  SEISMOLOCICAL  PROBLEMS* 

A Dynamic  Finite  Element  Method  (DFEM)  for  elastic  wave 
propagation  in  general  two-dimensional  heterogeneous  media 
was  developed  from  a causal,  variational  statement  of  Hamil- 
ton's Principle  and  applied- to  several  seismological  prob- 
lems. Judged  from  test  cases  where  the  results  were  compared 
against  Known  analytic  solutions,  the  method  was  able  to 
reproduce  temporal  and  phase  characteristics  of  Rayleigh 
waves  such  as  travel  time,  qroup  velocity  dispersion,  ellip- 
ticity,  and  phase  lag  to  within  5\  of  the  values  predicted  by 
theory.  The  Rayleigh  wave  amplitude,  however,  was  found  to 
be  25%  small  over  the  whole  spatial  waveform,  an  effect  which 
was  attributed  to  the  artificial  viscosity,  in  the  method 
necessary  to  remove  severe  numerical  noise.  Three  seismolog- 
ical  problems  of  interest  were  solved:  1)  Rayleigh  wave 
propagation  across  a continental  margin  structure,  2)  radia- 
tion from  an  extended  shear  source,  and  3)  a static  and 
dynamic  study  of  the  February  Q,  1971  San  Fernando  earth- 
quake. Conspicuous  among  the  many  results  were;  phase 

velocity  anisotropy  on  the  order  of  5%  in  the  continental 
margin  problem  which  aqreed  with  Alexander's  ( 1^63)  model 
experiment,  P and  S wave  radiation  spectra  from  the  extended 
shear  source  problem  which  fit  the  prediction  of  the  Brune 

(1970)  seismic  source  theory,  and  static  and  dynamic  models 

of  the  San  Fernando  earthquake  which  supported  the  results  of 

* 

Ph.D.  Thesis  of  D.  W.  McCowan 


previous  investigations  of  .near-field  data  for  this  event 
(e.g..  Bolt,  1972;  Alevine  and  Jordan,  1973)  . Based  on  these 
results,  conclusions  were  drawn  concerning  the  extension  of 
the  eethod  to  more  complicated  seisaological  problems. 
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ABSTRACT 


DETERMINATION  OF  SEISMIC  SOURCE  PARAMETERS  USING 
FAR- FI ELD  SURFACE  WAVE  SPECTRA* 

Far-field  surface  wave  amplitude  spectra  have  been 
studied  for  their  use  in  determining  seismic  source  parameters.  Using 
theoretical  spectra  generated  for  a double  couple  source  in  a plane 
layered  half-space,  the  minimum  residual  fit  to  the  observed  spectra 
was  obtained,  yielding  source  parameter  estimates. 

An  iterative  least-squares  regression  program  has  been 
developed  to  fit  observed  Rayleigh  and  Love  wave  amplitude  spectra. 
Iteration  is  on  the  dip  angle,  slip  angle,  strike  direction,  and  seismic 
moment.  The  initial  choice  of  these  parameters,  with  nominal  azi- 
muthal  coverage  (six  stations  over  180°  azimuth),  was  found  to  be  pri- 
marily limited  by  the  symmetry  of  the  theoretical  source  description. 
As  determined  from  tests  with  theoretical  spectra,  this  initial  choice 
can  be  as  much  as  30°  from  the  correct  value  of  the  angular  paramet- 
ers and  several  times  the  value  of  the  moment.  Intelligent  estimates 
of  these  initial  choices  can  be  obtained  using  either  bodywavc  source 
mechanism  estimates,  or  an  improved  version  of  Tsai  and  Aki's 
(1970a)  exhaustive  fitting  procedure  developed  in  this  study.  One  ad- 
vantage of  the  iterative  method  is  that,  for  the  same  parameter  set 
and  the  same  degree  of  precision  desired,  it  is  several  limes  faster 
computationally  than  exhaustive  fitting  procedures. 
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Aii  important  ex  tea  s io.n  of  both  the  iterative  and  ex- 
haustive techniques  that  has  been  developed  in  this  study  is  the  ampli- 
tude spectral  ratio  concept.  A 'master  event'  is  determined  for  a 
particular  geographical  region,  whose  source  mechanism  has  been 
established  with  a high  degree  of  confidence.  The  source  mechanism 
of  other  events  occurring  in  the  same  region  can  then  be  determined 
using  amplitude  spectral  ratios  with  the  master  event,  eliminating  the 
need  to  make  path  and  instrument  corrections. 

Complete  expressions  for  the  next  two  higher  multipoles 
(cjuadrapole  and  octapole)  were  generated  from  the  work  of  Harkrider 
(1970).  Their  importance  relative  to  that  of  a double  couple  was  found 
to  be  a function  of  frequency  and  the  physical  extent  of  the  rupture 
volume.  Alternatively,  an  extended  fault  can  be  modelled  using  multi- 
ple double  couples  displaced  in  both  space  and  time  which  simulate  a 
finite  rupture  velocity.  This  source  description  can  produce  non- 
symmetrical  surface  wave  radiation,  but  for  most  earthquakes  with 
magnitudes  less  than  6.  0 observed  at  teleseismic  distances,  a 
point  double  couple  source  description  is  adequate. 

Because  the  spectral  fitting  methods  depend  upon  the 
shape  and  level  of  the  amplitude  spectra,  theoretical  fundamental  and 
first  higher  mode  Rayleigh  and  Love  wave  double  couple  and  quadra- 
pole  amplitude  spectra  were  generated  for  many  source  configuration:;. 
The  most  important  observations  from  these  spectra  arc;:  (1)  for  both 


source  types  and  modes,  depth  variations  produce  the  largest  changes 
in  spectral  levels  and  shapes,  (2)  for  both  modes,  the  double  couple 
spectral  level  is  generally  an  order  of  magnitude  greater  than  that  of 
the  quadrapole,  (3) for  both  source  types,  the  largest  changes  in 
spectral  level  occur  in  the  10  to  20  seconds  period  range  for  the  fun- 
damental mode  and  9 to  12  seconds  period  range  for  the  first  higher 
mode,  and  (4)  the  most  prominent  spectral  holes  of  the  two  Rayleigh 
modes  occur  for  either  source  types  only  when  the  fault  is  vertical 
strike-slip. 

In  order  to  develop  procedures  for  obtaining  improved 
surface  wave  magnitude  estimates,  various  combinations  of  surface 
wave  spectra  for  the  fundamental  and  first  higher  mode  were  examin- 
ed, using  a double  couple  source,  to  assess  invariance  with  respect  to 
source  parameter  variations.  From  the  spectral  combinations,  the 
most  important  observations  are:  (1)  for  all  except  very  shallow 
sources  ( «•  5 km),  the  most  reliable  estimates  of  surface  wave  magni- 
tude ai  c obtained  from  azimuthal  averages  of  total  surface  wave 
2 2 

energy,  R + L.  , of  the  fundamental  mode,  and  (2)  for  shallow  sources 
it  is  best  to  use  an  azimuthal  average  of  fundamental  mode  Rayleigh 
waves  alone,  11  , in  the  period  range  of  29  to  30  seconds. 

Finally,  from  the  analysis  of  several  earthquakes  using 
both  spectral  fitting  procedures  and  spectral  ratios,  the  major  con- 
clusion.; are:  (1)  the  source  region  earth  model  has  little  . ffect  on  the 


depth  estimate  for  shallow  events  (<-10  km),  (2)  the  u/.imuthal  spread 
of  stations  should  ho  greater  than  90°  in  order  to  obtain  a valid  solu- 
tion, (3)  for  Eurasian  events  with  travel  path;,  along  the  Alpine- 
Himalavau  fold  system,  abnormally  low  group  velocities  (~0.  5km /sec 
low)  and  almost  twice  tile  normal  energy  attenuation  coefficient  were 
observed,  (4)  for  the  Bear  Valley  earthquake  of  22  June  1973,  the  far- 
field  seismic  moment  is  an  order  of  magnitude  greater  than  that  ob- 
tained from  near-field  observations,  and  (5)  that  the  use  of  spectral 
ratios  was  found  to  be  a valid  approach  in  the  analysis  of  two  pairs  of 
events. 
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The  Nature  and  Origin  of  Seismic  Codas  from  NTS  Explosions  Recorded  at  NORSAR 

On  the  basis  of  the  travel-time  tables  and  theoretical  studies  of  explosive 
sources,  short-period  seismograms  of  NTS  explosions  recorded  at  NORSAR,  an  epi- 
central  distance  of  approximately  73°,  should  show  relatively  simple  waveforms. 

The  first  arrival,  consisting  of  the  direct  P-wave  and  the  surface  reflection 
pP,  should  be  followed  18  seconds  later  by  PcP,  also  contaminated  by  its  surface 
reflection,  pPcP.  This  is  in  turn  followed  by  the  reflected  phases,  PP,  at  161 
seconds  after  P,  and  PPP  sane  263  seconds  after  P.  Each  of  these  arrivals  could 
reasonably  be  expected  to  consist  of  2 or  3 cycles  of  a predominantly  1 Hz  wave 
when  recorded  by  the  NORSAR  short  period  vertical  seismometers.  Figures  1 and 
2 show  typical  NORSAR  beams  from  explosions  detonated  in  two  distinct  areas  of 
NTS,  Pahute  Mesa  and  Yucca  Flats  respectively.  The  pertinent  epicenter  data  for 
these  events  are  given  in  Table  1.  Clearly  these  seismograms  are  far  from  simple, 
particularly  those  from  the  Pahute  Mesa  events  (Figure  1).  This  observation 
prompts  the  obvious  question  as  to  why.  In  the  discussion  that  follows  we  will 
attempt  to  elucidate  the  nature  of  the  seismic  coda  and  determine  the  source  of 
its  generation. 

One  of  the  valuable  properties  of  an  array  is  its  ability  to  resolve  the 
speed  and  propagation  direction  of  an  arrival.  The  array  can  then  be  steered  in 
this  direction  and  beams  formed  so  as  to  minimize  the  contamination  of  the  desired 
signal  by  unwanted  arrivals.  Not  only  can  the  expected  arrivals  be  separated  on 
the  basis  of  arrival-time  at  a single  station,  they  can  be  separated  on  the  basis 
of  propagation  velocity  across  the  array.  A systematic  approach  is  to  assume  that 
all  signals  originating  in  the  source  region  arrive  at  the  receiver  along  the  same 
azimuth  as  the  direct  P waves,  and  to  form  beams  corresponding  to  the  velocity 


range  of  the  expected  late  arrivals.  The  relative  power  in  each  beam  for  a given 


window  is  then  displayed  as  a function  of  beaming  velocity  and  arrival  time. 
Contours  of  equal  power  then  reveal  coherent  arrivals  as  peaks  in  velocity  vs 
time  space. 

The  lower  part  of  Figure  3 shows  such  a plot  for  the  explosion  STILTON, 
which  was  detonated  in  Pahute  Mesa,  slightly  to  the  NW  of  the  Silent  Canyon 
Caldera  (see  Table  1);  the  contour  interval  is  3 db.  The  first  arrival  has  an 
apparent  velocity  of  18.0  km/sec,  which  is  somewhat  lower  than  the  18.8  km/sec 
predicted  by  the  Herrin  travel-time  tables,  possibly  due  to  dipping  interfaces 


within  the  crust  beneath  the  array.  However,  the  important  feature  of  this 
plot  is  that  the  conspicuous  arrival  18  seconds  after  the  initial  P-wave  onset, 
also  has  an  apparent  velocity  of  18.0  km/sec.  So,  although  the  differential 
travel-time  is  correct  for  PcP,  the  velocity  is  much  slower  than  the  PcP  pre- 
dicted value  of  25.8  km/sec.  In  fact,  all  the  arrivals  in  the  90  seconds  of  the 
beam  following  the  initial  P-wave  onset  have  an  apparent  velocity  of  18.0  km/sec. 
This  is  amply  demonstrated  by  normalizing  the  plot  by  the  maximum  in  each  time 
window,  as  is  shown  in  the  upper  part  of  Figure  3,  which  is  contoured  at  1 db 
intervals.  A search  for  arrivals  was  conducted  by  steering  the  array  toward  + 

20  degrees  from  the  P-wave  back-azimuth.  No  significant  off-azimuth  arrivals  were 


detected. 

Although  the  velocity  time  plots  are  themselves  a qualitative  measure  of  the 
coherency  of  signals  across  the  array,  a more  quantitive  measure  of  coherency  is 
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M is  the  length  of  the  time  window  over  which  S is  computed 
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Figure  4 shows  the  coherency  in  the  best-beam  from  STILTON  as  function  of  time; 
the  relative  power  in  the  beam  as  a function  of  time;  and  the  beam  itself.  They 
confirm  that  the  arrival  18  seconds  after  the  P-wave  has  the  same  coherency  as 
the  initial  P-wave  at  the  beaming  velocity  of  18  km/sec,  and  is  only  7 db  down 
in  power.  Thus,  this  arrival  cannot  be  PcP.  The  plots  show  that  further  dis- 
crete arrivals,  also  coherent  at  the  P-wave  velocity,  occur  up  to  60  seconds 
after  the  initial  P-wave  arrival. 

Beam  steering,  coherency  and  relative  power  plots  were  produced  for  each 
of  the  8 Pahute  Mesa  events  shown  in  Figure  1.  These  tended  to  confirm  the  results 
from  the  analysis  of  STILTON.  No  evidence  of  PcP,  or  other  phases  arriving  with 
velocities  and  propagation  directions  significantly  different  from  the  initial 
P-wave  were  found.  The  measured  P-wave  velocities  themselves  were  identical  to 
within  the  experimental  error.  The  distribution  in  time  of  the  coda  arrivals, 
their  relative  power,  and  their  coherencies,  all  varied  from  event  to  event.  As 
an  example.  Figure  5 shows  the  coherency  and  relative  power  from  KASSERI.  This 
event  clearly  has  many  more  arrivals  that  STILTON  and  the  coda  remains  both 
coherent  and  more  energetic  over  a longer  interval.  No  arrivals  were  found  with 
velocities  significantly  different  than  that  of  the  initial  P-wave. 
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To  examine  the  coherency  of  the  codas  among  events  at  Pahute  Mesa,  the  8 
NORSAR  beams  were  stacked.  Each  beam  was  normalized  to  unit  amplitude,  and  time 
aligned  on  the  front  trough  of  the  P-wave  signal.  The  lower  trace  in  Figure  6 
shows  the  result.  The  maximum  amplitude  in  the  P-wave  is  0.89,  so  that  the  signal 
loss  due  to  dissimilarities  in  waveform  (see  Figure  1)  is  not  severe.  In  this 
case  of  coda,  however,  the  amplitudes  are  reduced  by  approximately  a factor  of  3 
throughout  the  time  interval,  and  thus  are  mostly  incoherent  from  event  to  event. 

A similar  set  of  analyses  were  performed  on  the  8 Yucca  Flat  events  (Figure  2). 
The  direction  and  speed  of  propagation  (apparent  velocity)  of  this  initial  P-wave 
from  these  events  was  essentially  the  same  as  for  the  Pahute  Mesa  events.  The  coda 
also  arrived  with  the  same  apparent  velocity.  However,  as  Figure  7 shows,  the 
relative  power  and  coherency  of  the  coda  arrivals  are  much  lower.  Despite  the 
greater  degree  of  similarity  of  the  initial  P-waves,  there  is  still  significant 
variation  in  time  and  amplitude  of  distinct  coda  arrivals  from  event  to  event. 

This  is  demonstrated  by  the  upper  trace  in  Figure  6 which  shows  the  phased  sum  of 
the  beams.  The  amplitude  of  the  P-wave  is  0.97  in  this  case,  whereas  coda  is 
reduced  by  a factor  of  approximately  2.5,  again  suggesting  that  the  coda  arrivals 
are  uncorrelated  from  event  to  event.  There  was  no  indication  of  significant 
arrivals  with  velocities  other  than  that  of  the  initial  P-wave,  in  particular  PcP. 

Articles  in  the  literature  on  the  generation  of  seismic  codas  are  legion. 
Discounting  contributions  by  other  phases  such  as  pP,  sP,  PcP  etc.,  published  mech- 
anisms fall  into  the  following  categories: 

(i)  multiply  reflected  phases  originating  in  the  crust  and  upper  mantle 
under  the  source  and/or  the  receiver 

(ii)  reflections  and  refractions  of  P-waves  taking  place  out  of  the  dia- 
metrical plane 
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(iii)  P-wave  scattering  of  the  Chernov  type  at  the  source  and/or  receiver 

(iv)  mode  conversions,  namely  P-to-Rayleigh  waves  at  the  receiver  and 
Rayleigh  to  P-wave  at  the  source 

(v)  tectonic  release  including  aftershocks  and  other  source  related 
phenomena. 

Of  these  mechanisms,  the  nature  of  the  codas  themselves  as  discussed  elim- 
inates all  possibilities  with  the  exception  of  those  which  generate  coda  in  the 
immediate  proximity  of  the  source.  It  is  well  known  that  explosions  at  NTS, 
particularly  those  at  Pahute  Mesa,  generate  considerable  SH  wave  radiation  in 
some  instances.  If  coda  generation  is  due  to  tectonic  release,  those  events  with 
high  coda  levels  should  also  be  correlated  with  events  with  large  SH  excitation. 

Table  3 shows  the  Love  to  Rayleigh  wave  amplitude  ratios  from  seismograms  recorded 
at  5 SDOS  stations  within  the  continental  United  States,  for  the  8 Pahute  Mesa 
events  and  6 of  the  Yucca  Flat  events.  There  is  no  obvious  correlation  between 
events  having  high  coda  levels  and  events  with  large  L/R  ratios.  Furthermore,  the 
Pahute  Mesa  events  with  coda/P-wave  ratios  of  0.5  would  require  aftershocks  of 
approximately  magnitude  6 to  account  for  the  observed  coda  levels.  Although  we 
have  not  examined  the  very  near-in  records  for  these  events,  there  is  no  evidence 
in  the  literature  to  suggest  that  such  events  are  triggered  18  to  20  seconds  after 
similar  large  explosions. 

Therefore,  we  are  left  with  mechanisms  (i),  (iii)  and  (iv)  to  consider.  Further- 
more, reflections  in  a plane-layered  structure  can  also  be  eliminated,  because  such 
reflections  would  be  well-correlated,  though  possibly  time-shifted,  from  event  to 
event  within  the  two  source  regions.  We  are  inclined  to  discount  the  Chernov  mech- 
anism, because,  strictly  speaking,  this  theory  is  only  applicable  to  small  random 
variations  in  elastic  parameters.  Because  of  the  size  of  the  coda  arrivals,  a second- 
order  or  strong-scattering  theory  wuld  be  required.  Thus  we  are  left  multiple 


6 

reflections  from  curved  boundaries  which  lead  to  the  formation  of  caustics  of  the 
type  discussed  by  Hong  and  Helmberger  (in  press),  and  Rayleigh  to  P-wave  conversions 
in  the  immediate  vicinity  of  the  source. 

In  an  effort  to  further  resolve  this  problem,  seismograms  of  teleseismic 
events  recorded  by  SDCS  stations  located  within  NTS  are  being  examined.  Preliminary 
results  show  that  a Russian  event  at  Novaya  Zemlya  generated  unusually  complex 
seismograms  at  NT-NV,  which  is  within  Pahute  Mesa,  a much  simpler  one  at  0B2NV, 
located  in  the  nearby  Climax  Stock.  Filtering  of  these  and  other  seismograms  from 
deep  events  will  be  used  to  determine  whether  these  codas  are  primarily  body  or 
surface  waves.  In  any  case  these  data  provide  strong  additional  support  for  the 
conclusion  that  scattering  near  the  source  is  the  primary  mechanism  for  coda  gen- 
eration associated  with  explosions  at  NTS.  Structural  complexities  beneath  Pahute 
Mesa  are  especially  effective  compared  to  other  test  sites  in  generating  large  and 
variable  coda  levels. 


Epicenter  Data  For  Pahute  Mesa  Fvents 
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Table  3 

PAHUTE  MESA  L/R  RATIOS 


Event/Sta. 

RK-ON 

CPS0 

WH2YK 

FN-WV 

HN-ME 



CODA 

C/P 

Kasser i 

0.39 

O.Sb 

0.87  . 

1.12 

6.74 

6.4 

H 

0.5 

Muenster 

0.82 

1.31 

0.99 

4.71 

6.2 

H 

0.4 

Pool 

0.90 

1.06 

0.90 

3.86 

6.1 

L 

0.25 

Cheshire 

0.43 

0.48 

0.73 

0.78 

1.85 

6.0 

L 

0.2 

Fontina 

0.42 

0.92 

0.61 

1.93 

6.3 

H 

0.3 

Colby 

0.16 

0.28 

0.18 

0.26 

6 . 3 

II 

0.5 

Inlet 

0.87 

1.02 

0.63 

1.28 

1.56 

6.0 

L 

0.2 

Estuary 

0.30 

0.62 

0.72 

0.70 

1.87 

6.0 

H 

0.45 

This  Table  illustrates  variability  of  tectonic  release  among  events  in  Pahute  Mesa 
and  the  lack  of  correlation  between  large  tectonic  release  and  teleseismic  coda 
levels  as  observed  at  NORSAR.  H denotes  high  coda  levels  and  L denotes  low  coda  levels. 


Time  (Seconds) 
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elative  power  as  a function  of  time  and  beaming  velocity  for  STILTON.  Power  plot  is  contoun 
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nly  -1  and  -3  db  contours  are  shown.  Figure  4c  shows  the  beamed  seismogram  at  NORSAR  for  ST 
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INTRODUCTION 

Conventional  body  wave  fault  plane  solutions  for  a seismic 
event  inherently  contain  considerable  scatter.  A wide  azimuthal 
distribution  of  body  wave  first  arrivals  is  generally  necessary 
in  order  to  minimize  this  scatter.  Fault  plane  solutions  for 
events  of  the  August,  1975  Oroville,  California  sequence  have 
been  determined  by  several  groups  of  investigators.  These  solutions 
are  both  numerous  and  overlapping. 

The  main  event  (ML  = 5*7)  of  the  Oroville,  California  se- 
quence occurred  in  a region  of  historically  low  seismicity 
approximately  10  kms . south  of  the  4 billion  cubic  meter  man- 
made Lake  Oroville.  The  sequence  began  seven  years  after  the 
onset  of  filling  but  only  one  month  after  its  most  rapid  filling. 

The  largest  event  prior  to  this  sequence  was  an  event  in  1940, 
also  a magnitude  5*7.  about  50  kms.  north  of  the  town  of  Oroville. 
Twenty-one  foreshocks  (M{(  > 1.6)  preceded  the  main  event  beginning 
on  June  28,  1975*  There  were  approximately  300  aftershocks 
recorded  by  local  stations  and  arrays  from  an  area  14  by  10 
kms.  southeast  of  Oroville  through  August  11,  1975* 

The  purpose  of  this  investigation  was  to  more  clearly 
define  the  focal  mechanisms  of  the  layer  events  of  thi3  sequence. 

The  apDroach  utilized  was  the  synthesis  of  long  period  Love  and 
Rayleigh  waves  whose  time  domain  signatures  and  frequency  domain 
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spectra  were  then  compared  with  the  regionally  observed 
signatures  and  spectra  recorded  at  the  Seismic  Research 
Observatory  (hereafter  SRO)  in  Albuquerque,  New  Mexico.  The 
focal  parameters,  namely  source  depth,  strike  azimuth,  fault 
plane  dip  and  slip  were  perturbed  in  order  to  obtain  the 
best  correlation  between  the  synthetic  and  observed  surface 
wave  signatures  and  spectra. 

Two  fault  plane  solutions  were  determined  which  would 
synthesize  the  observed  signatures  and  spectra  fairly  well  as 
recorded  at  the  SRO  station  in  Albuquerque.  One  solution 
corresponded  to  the  sequence  of  foreshocks,  main  event,  and 
some  aftershocks,  whereas,  the  second  solution, which  had  a 
distinctly  different  source  depth  and  orientation,corresponded 
to  the  aftershock  on  August  6,  1975* 

METHODS 

A suite  of  ten  events  (ML>  4*3)  which  were  digitally 
recorded  at  the  SRO  station  in  Albuquerque  was  obtained  for 
the  Oroville,  California  sequence.  This  suite  consisted  of 
foreshocks  to  the  main  event  on  August  1,  1975.  the  main  event, 
and  the  aftershocks  through  August  6,  1975*  Locations  for 
some  of  these  events  are  shown  in  the  map  of  the  Oroville 
area  of  Figure  1.  These  long  period  digital  recordings  were 


sampled  at  one  sample  per  second.  These  digital  recordings 
were  amenable  to  rather  easy  computer  manipulations  such  as 
time  scale  changes,  gain  changes,  rotation  and  Fourier  analysis. 


Figures  2 through  4 are  the  vertical,  radial  and  transverse 
components,  respectively,  of  four  of  these  events  recorded  at 
Albuquerque.  The  scale  factors  are  indicated  by  each  trace. 
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. 


Note  that  the  5*8  mainshock  is  clipped. 

The  synthetic  seismograms  were  generated  assuming  a 
plane  parallel  layered  source  structure.  A three  step 
computational  scheme  calculated  the  seismograms  containing  root 
periods  ranging  from  5 to  75  seconds.  The  resulting  vertical, 
radial  and  transverse  component  seismograms  were  then  compared 
visually  to  the  rotated  seismograms  recorded  at  Albuquerque. 

The  comparison  was  made  by  overlaying  the  observed  signatures 
with  the  synthetic  signatures,  noting  the  goodness-of-f it 
of  the  relative  amplitudes  and  duration  of  the  individual  wavelets 
comprising  the  seismograms.  In  addition,  the  Love  and  Rayleigh 
wave  spectral  ratios  were  compared  as  a function  of  frequency 
for  the  synthetic  and  observed  events. 

PROCEDURE  FOR  SYNTHESIS 

The  Fortran  IV  computer  programs  employed  in  the  synthesis 
of  surface  waves  propogated  from  Oroville,  California  to 
Albuquerque,  New  Mexico  are  the  implementation  of  the  theory 
developed  by  Harkrider  (1964),  Ben-Menahem  and  Harkrider  (1964) 
and  Harkrider  (1970).  The  implementation  required  two  runs  of 
a series  of  three  programs  for  the  synthesis  of  Rayleigh  and  Love 
wave  seismograms.  A description  of  the  input  card  formats  can 
be  found  in  the  appendix  of  Alexander,  et  al.  (1973). 
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An  important  assumption  which  was  required  to  be  made  in 
this  study  wa3  the  specification  of  a suitable  source  structure 
for  the  Oroville  area.  The  structure  consisted  of  model  35  CM2 
(Alexander,  1963)  with  the  upper  33  kms . replaced  by  the 
structure  determined  by  refraction  profiles  (Eaton,  1966). 

This  upper  section  reflects  the  local  geologic  environment  of 
valley  sediments  overlying  Jurassic-Triassic  metavolcanics 
(Ryall  and  Van  Wormer  1975)  in  the  Oroville  region.  The  structure 
used  in  this  study  is  shown  in  Figure  5 and  tabulated  in  Table  1. 

The  dispersion  programs  were  run  first.  They  calculated 
the  depth  independent  quantities  such  as  phase  velocities, 
group  velocities  and  the  amplitude  response  as  a function  of  root 
period  for  Love  and  Rayleigh  waves.  The  surface  ellipticity 
v,ii  also  calculated  in  the  case  of  Rayleigh  waves.  The  required 
input  at  this  stage  was  the  local  source  structure  as  given 
above  as  well  as  the  specification  of  the  periods  at  which  roots 
were  to  be  found.  The  periods  ranged  from  5 to  75  seconds  in 
this  study. 

The  output  from  the  dispersion  programs  was  used  as  the 
input  to  the  punch  programs.  The  punch  programs  calculated  the 
displacement-stress  quantities  at  specified  depths  in  the  source 
structure.  The  focal  depths  investigated  in  this  study  were  2, 

3,  5,  7,  8,  9,  10  and  12  kms.  which  were  consistent  with  the 
range  of  focal  depths  determined  for  the  Oroville  sequence  by 
various  investigators  (see  Table  2). 

Outputs  from  both  the  dispersion  and  punch  programs  were 
then  used  as  inputs  to  the  step  which  computes  the  Rayleigh 
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and  Love  wave  seismograms  in  a plane  parallel  layered  structure. 

The  synthetic  seismograms  that  were  computed  were  corrected  for 
the  instrument  response  of  the  SRO  station  in  Albuquerque  as 
shown  in  Figure  6.  The  attenuation  coefficients  employed  were 
determined  by  Hermann  and  Mitchell  (1975)  for  the  interior 
of  North  America.  These  coefficient  values  contained  considerable 
scatter  and  were  therefore  used  only  as  a first  approximation 
of  the  attenuation  properties  of  the  western  United  States. 

The  values  used  in  this  study  are  tabulated  in  Table  3» 

The  additional  input  parameters  utilized  in  the  final  step 
of  the  synthesis  were  those  corresponding  to  the  fault  plane 
orientation,  source  type,  source  time  function  and  cpicentral 
distance.  The  fault  plane  orientations  were  specified  following 
the  convention  of  Ben-Menahcm  and  Harkrider  (1964).  Synthetic 
Love  and  Rayleigh  waves  were  computed  utilizing  a source  function 
specified  by  either  a point  source  or  a double  couple  source. 

For  the  case  of  the  double  couple,  the  source  time  function  was 
employed  corresponding  to  a magnitude  5 event  determined  from 
Aki  (1967). 

EXPERIMENTS  AND  RESULTS 

Numerous  Rayleigh  and  Love  wave  seismograms  modeling  the 
Oroville,  California  events  were  synthesized  using  the  aforementioned 
procedure.  Fublished  fault  plane  solutions  of  foreshocks,  the  main 
event  and  aftershocks  listed  in  Table  2 served  as  a starting  point 
for  the  solutions  for  the  synthetic  wave  forms.  The  trial  fault 
plane  solutions  employed  in  the  synthesis  were  then  perturbed  to 


cover  the  ranges  of  the  published  values. 

The  vertical,  radial  and  transverse  component  synthetic 
seismograms  were  visually  compared  with  a rotated  event 
recorded  at  Albuquerque,  New  Mexico,  both  plotted  with  the  same 
time  scale (Figure  7)*  This  event  was  a magnitude  4.7  (BRK) 
foreshock  four  hours  prior  to  the  August  1,  1975  magnitude  5.7 
(3RK)  Oroville  sequence  main  event.  The  main  event  was  not 
used  for  comparison  in  this  study  for  two  reasons.  The  surface 
wave  coda  of  the  main  event  was  contaminated  by  the  coda  of  a 
foreshock  eight  seconds  prior  to  the  main  event.  Tn  addition, 
the  wave  form  as  recorded  at  Albuquerque  was  clipped  (see  Figures 
2 through  4).  This  comparison  of  the  synthetic  and  observed 
vertical,  radial  and  transverse  components  permitted  the  rejection 
of  many  solutions  as  viable  mechanisms  for  the  observed  wave  forms. 
It  was  found  that  the  depth  of  focus  was  not  as  sensitive  in  these 
eyeball  matches  of  synthetic  versus  real  waveforms  as  were  the 
strike,  dip  and  slip  of  the  particular  fault  plane  solution 
under  scrutiny. 

The  depth  of  focus  was  localized  by  comparing  the  Love  to 
Rayleigh  wave  spectral  ratios  for  the  synthetic  and  observed 
data  as  seen  in  Figure  8.  The  spectra  of  the  observed  data  was 
obtained  by  employing  the  usual  Fourier  techniques.  The  peak 
in  the  synthetic  Love  to  Rayleigh  wave  ratio  was  observed  to  shift 
to  lower  frequencies  a3  the  depth  of  focus  was  increased. 

The  best  fitting  synthetic  seismogram  employed  a double 
couple  extended  source  function.  The  fault  plane  parameters 
which  yielded  this  best  fit  with  the  foreshock  were  azimuth,  N 10PWj 
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dip,  60'Wj  slip,  90".  The  slip  angle  of  90"  is  consistent 
with  the  angle  listed  in  the  figures,  270  . This  angle  is 
specified  in  the  convention  of  Ben-Menahem  and  Harkrider  (1964) 
and  represents  a normal  dip  slip  mechanism.  In  this  case  the 
Great  Valley  block  has  moved  down  with  respect  to  the  Sierra 
Nevada  block.  The  depth  determined  by  spectral  ratios  was 
5 kilometers.  The  epicentral  distance  from  Oroville  to  Albuquerque 
used  in  this  analysis  was  1426  kms.  The  Rayleigh  and  Love  waves 
computed  with  this  fault  plane  solution  also  compared  favorably 
with  several  other  aftershocks  of  the  Oroville  sequence.  It 
was  therefore  concluded  that  the  main  event  could  also  be 
characterized  by  this  fault  plane  solution. 

Further  examination  of  the  aftershock  data  for  the  Oroville 
events  revealed  some  marked  changes  in  the  recorded  waveforms 
(see  Figures  2 through  4).  The  changes  in  the  observed  waveforms 
for  the  August  6,  1975  aftershock  were  adequately  modeled 
utilizing  a different  fault  plane  solution  than  that  determined 
above.  A sot  of  approximate  fits  was  chosen  from  the  original 
suite  of  synthetic  seismograms.  Again,  the  depth  was  localized 
using  the  Love  to  Rayleigh  spectral  ratios.  The  focal  parameters 
which  best  modeled  the  August  6 aftershock  were  azimuth,  N 0®  Wj 
dip,  40°Wj  slip,  90^i  and  depth  10  kms.  These  results  can  be 
seen  in  Figures  9 and  10. 

DISCUSSION 


The  results  of  this  study  indicate  that  it  is  possible 
to  determine  the  source  mechanisms  of  regional-teleseismic  events 
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utilizing  broad  band  data  recorded  at  only  a few  high  quality 
stations.  The  Seismic  Research  Observatory  network,  of  which 
there  are  ten  installations  in  current  operation,  supply  such 
high  quality  data.  In  this  study  Love  and  Rayleigh  wave  data 
from  only  one  station  were  used. 

The  fault  plane  solutions  determined  in  this  study  were 
found  to  be  in  good  agreement  with  the  range  of  solutions 
published  previously.  The  change  in  focal  parameters  between  the 
foreshock  and  aftershock  of  August  6 could  be  associated  with 
the  development  of  a non-uniform  local  stress  distribution 
following  the  mainshock.  The  method  employed  in  this  study  was 
found  to  be  quite  sensitive  to  small  changes  in  fault  plane 
parameters  as  seen  in  the  two  additional  waveforms  in  Figure  9. 
Therefore,  this  approach  may  be  of  significant  help  in  determining 
fault  plane  solutions  of  regional  events  when  station  coverage 
is  poor,  provided  that  an  appropriate  source  structure  can  be 
provided . 

In  conclusion,  it  must  be  emphasized  that  the  proper 
propagation  structure  should  be  employed  in  the  synthesis  of 
surface  waves.  This  study  revealed  the  rather  fortuitous 
result  that  the  propagation  path  structure  was  identical  to  the 
source  structure  given  in  Figure  5*  Additional  theoretical 
seismograms  shown  in  Figures  11  and  12  which  did  not  match  the 
observed  waveforms  were  synthesized  using  the  focal  parameters 
determined  above.  The  source  structure  remained  the  same  as 
used  previously,  but  the  propagation  path  was  that  of  35  CM2 
(Alexander,  1963).  It  is  quite  evident  that  this  change  does  not 
produce  satisfactory  results. 
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Table  1 


Oroville  Structure 


Thickness,  kms. 


Alpha  km/sec.  Beta  km/sec.  Rho  gm/cm ^ 


2.0 

3.0 

5.0 
5.0 

5.0 
5-0 
R.O 

8.0 
8.0 
7.5 
7.5 

10.0 

10.0 

10.0 

10.0 

10.0 

10.0 

20.0 

20.0 

20.0 

20.0 

100.0 


3-75 

5.90 

6.80 

6.80 

6.80 

6.80 

6.80 

8.00 

8.00 

8.09 

8.09 

7.96 

7.96 
7.84 
7.84 
7.82 
7.82 
7.86 
7-9  2 

7.97 
8.03 
8.51 


1.72 

3.40 
3.88 
3.88 
3-88 
3.88 
3.88 

4.40 

4.40 
4.65 

4.65 
4.56 
4.56 
4.48 
4.48 
4.43 
4.43 

4.41 

4.39 

4.40 

4.41 

4.66 


2.20 

2.80 

3.03 

3.03 

3.03 

3.03 

3.03 

3.26 

3.26 

3-37 

307 

3.42 

3.42 

3.43 

3.43 

3.44 

3.44 

3.45 
3-46 
3-48 
3.52 
3-58 


Fault  Plane  Solutions  of  the  Oroville  Sequence 
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Table  3 

Attenuation  Coefficients  used  in  the  Synthesis  * 


Period 

Rayleigh 

xl0-4/km 

Love 

xl0“4/km 

75 

2.10 

1.20 

70 

2.00 

1.15 

*5 

1.90 

1.05 

60 

1.75 

0.95 

55 

1.60 

0.90 

50 

1.5  0 

0.80 

45 

1-35 

0.70 

40 

1.3  0 

0.60 

35 

1.3  0 

0.55 

30 

1.25 

0.50 

25 

1.25 

0.40 

20 

I.25 

O.30 

19 

I.30 

0.25 

18 

1-35 

0.25 

17 

1.40 

0.20 

16 

1.45 

0.25 

15 

1.50 

0.25 

14 

1.60 

0.35 

13 

1.80 

0.45 

12 

1.90 

0.65 

11 

2.10 

1.20 

10 

2.30 

2.30 

9 

2.70 

3.50 

0 

3.30 

4.70 

7 

5.00 

5.90 

6 

6.10 

7.00 

5 

7.20 

7.60 

* From  Hermann  and  Mitchell  (1975) 
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Figure  Captions 


Figure  1 


Figure  2 


Figure  3 


Figure  4 


Location  map  of  the  Oroville,  California  region 
indicating  some  of  the  epicenters  of  the  sequence. 

From  Morrison  et  al.  (19?6). 

Long-period  vertical  component  signals  recorded  at 
the  Albuquerque,  New  Mexico  SRO  station  (AMNO)  for 
4 events  from  Oroville,  California  in  August,  1975* 
Relative  signal  strength  is  indicated  by  scale 
factors  Z max.  The  mainshock  (bottom  trace)  is  clipped. 

Long-period  radial  component  signals  recorded  at 
the  Albuquerque,  N.  M.  SRO  station  (AMNO)  for  4 events 
from  Oroville,  California  in  August,  1975*  Relative 
signal  strength  is  indicated  by  scale  factors  R max. 

The  mainshock  (bottom  trace)  is  clipped. 

Long-period  transverse  component  signals  recorded  at 
the  Aubuquerque,  N.  M.  SRO  station  (AMNO)  for  4 events 
from  Oroville,  California  in  August,  1975*  Relative 
signal  strength  is  indicated  by  scale  factors  T max. 

The  mainshock  (bottom  trace)  is  clipped. 


Figure  5 

Figure  6 
Figure  7 

Figure  8 

Figure  9 


Plot  of  the  Oroville  structure  employed  for  the 
synthetic  seismograms  of  this  study.  Values  are 
tabulated  in  Table  1. 

Amplitude  response  of  the  Seismic  Research  Observatory 
station  in  Albuquerque,  New  Mexico. 

Comparison  of  the  observed  (SRO)  and  theoretical 
signals  at  Albuquerque,  N.  M.  for  the  August  1,  1975 
Oroville  foreshock.  Source  parameters  for  the 
theoretical  waveforms  are  given  in  Figure  8. 

Comparison  of  observed  and  theoretical  Love  to  Rayleigh 
spectral  ratio  (L/R)  for  the  signal  shown  in  Figure  7. 
The  legend  gives  the  source  parameters  (depth,  strike, 
dip  and  slip)  for  the  theoretical  waveform  synthesis. 

Comparison  of  the  observed  (SRO)  and  theoretical 
signals  at  Albuquerque,  N.  M.  for  the  August  6,  1975 
Oroville  aftershock.  Source  parameters  for  the 
theoretical  waveforms  shown  are  given  in  Figure  10. 


Figure  10  Comparison  of  observed  and  theoretical  Love  to  Rayleigh 
spectral  ratio  (L/R)  for  the  signal  shown  in  Figure  9. 
The  legend  gives  the  source  parameters  (depth,  strike, 
dip  and  slip)  for  the  theoretical  waveform  synthesis. 
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Figure  11 


Figure  12 
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Synthetic  seismogram  with  same  focal  parameters 
and  source  structure  as  Figure  ?•  Propogation  path 
structure  is  35  CM2.  Compare  with  Figure  7. 


Synthetic  seismogram  with  same  focal  parameters 
and  source  structure  as  Figure  9,  Propogation  path 
structure  is  1*5  CM2.  Compare  with  Figure  9. 
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Abstract 


. 1WP  SJl 


We  derive  a shear  wave  crust  and  upper  mantle  structure  for  the 
southern  part  of  Novaya  Zeralya  by  an  application  of  the  two  event, 
single  station  method  of  Rayleigh  wave  phase  velocity  dispersion 
analysis.  This  method  provides  a means  of  isolating  the  surface  wave 
dispersion  characteristics  of  a remote  source  region  using  only 
teleseismic  recordings.  The  observed  phase  velocity  data  are  then  sys- 
tematically inverted  to  obtain  a best-fitting  model.  Our  preferred 
model  has  a 45  km  thick  crust  with  no  shear  wave  low  velocity  zone  in 
the  upper  mantle.  It  is  similar  to  published  structures  for  the  south- 
ern Ural  mountains  and  is  therefore  compatible  with  the  premise  that 
Novaya  Zemlya  is  a northern  extension  of  the  Urals. 
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Introduction 


A sequence  of  large,  presumed,  explosions  detonated  at  Novay a 
Zemlya  between  1972  and  1975  generated  unusually  good  Rayleigh  wave 
seismograms.  A set  of  five  recordings  of  these  explosions  made  at  the 
Alaskan  Long  Period  Array  (ALPA)  was  used  in  an  attempt  to  infer  the 
crust  and  upper  mantle  structure  of  the  source  region.  It  appears 
from  the  epicenter  information  given  in  Table  1 that  there  are  two 
distinct  test  sites  on  Novaya  Zemlya.  The  furthest  from  ALPA  is  at 
the  southern  end  of  the  lower  island  (see  Figure  3).  The  nearer  test 
site  is  located  at  the  fiord  separating  the  lower  island  from  the  rest 
of  Novaya  Zemlya.  The  average  distance  between  the  two  test  sites  is 
285  kra. 

Fortuitously,  these  two  test  sites  lie  nearly  on  a great  circle 
path  to  ALPA  (shown  in  Figure  2).  A large  portion  of  that  path 
traverses  the  Arctic  ocean  basin  neatly  passing  its  way  between  the 
North  American  and  Eurasian  continental  shelves.  Since  the  propagation 
path  to  ALPA  traverses  relatively  simple  structures,  we  can  expect  the 
effects  of  anomalous  propagation  between  the  nearer  test  site  and  AL- 
PA, that  is  along  the  common  portion  of  the  path,  to  be  minimal. 

The  combination  of  the  above  factors  affords  the  unique  opportun- 
ity to  Isolate  the  phase  velocity  dispersion  due  to  propagation  in  the 
source  region  itself.  The  method,  due  to  Alexander  (1969)  and  used  by 
Taylor  (1972)  to  study  the  mid-Atlantic  ridge,  is  a modification  of 
the  usual  two  station,  single  event  procedure  (e.g.,  Glover  & Alex- 
ander, 1969).  However  here,  instead  of  finding  the  phase  velocity 
dispersion  between  the  two  stations,  the  combination  of  me  station 
and  two  events  facilitates  finding  the  phase  velocity  dispersion  in 


the  source  region. 


Method 

We  first  computed  the  array-averaged  cross  power  spectrum  for 
different  north-south  pairs  of  events  at  Novaya  Zemlya  using  a group 
velocity  window  of  2.6  to  4.0  km/sec.  Symbolically  this  is. 


R(t)-“f  IX  (t)Y  (t+r) 

" J=1  T 0 J 

P(»)  = | R(t)^iutdt  (1) 

Here  X and  are  the  seismograms  from  the  event  pair  recorded  at  the 
ALPA  sensor,  R is  the  array-averaged  cross  correlation  function, 
and  P is  the  resulting  cross  power  spectrum.  This  can  then  be  ex- 
pressed as  an  amplitude  and  phase: 

P(w)  -|P(a))|ei<5(w)  (2) 

If  the  original  seismograms  contain  common  instrument  and  source 
phases,  these  cancel  when  computing  the  cross  power  spectrum.  |P(ui)f 
then  becomes  the  average  cross  power  and  <S(u$  is  related  to  the  phase 
velocity  in  the  source  region,  C(u>),  by: 

tuD  ,, 

«">  • Z0-  ±2"n  <3> 

where  D ^ is  the  projection  of  the  event  separation  on  the  great  cir- 
cle path  shown  in  Figure  1.  Our  average  event  separation,  285  km,  al- 
lowed us  to  pick  the  arbitrary  integer  easily  by  trial  and  error.  An 
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earth  structure  model  is  then  obtained  by  matching  its  phase  velocity 
dispersion  to  that  obtained  by  the  data  analysis* 


Implicit  in  this  analysis  are  the  following  four  assumptions: 

1*  The  source  region  can  be  characterized  by  an  average  plane 

layered  structure 

2.  The  Rayleigh  waves  from  all  events  traverse  a common  path 

from  the  nearer  (northern)  events  to  the  ALPA  array 

3.  Interference  due  to  multipathing  can  be  neglected 

4.  The  source  phase  for  the  two  events  is  identical 

Data  Analysis 

For  this  study  we  had  five  events,  three  northern  and  two  south- 
ern. The  pertinent  epicenter  information  is  given  in  Table  1 and  ALPA 
center  sensor  seismograms  are  shown  in  Figure  4.  As  can  be  seen, 
these  are  all  large  explosions  which  excited  sufficiently  strong  Ray- 
leigh waves,  despite  their  non-slip  source  mechanisms,  so  that  low 
signal-to-noise  ratios  were  not  a problem.  In  order  to  rtaove  loca- 
tion bias,  all  five  events  were  relocated  using  P-vaves  observed  at  a 
common  set  of  35  stations.  After  constraining  the  depth  of  each  event 

to  zero,  the  areas  of  the  resulting  epicenter  confidence  error  el- 

2 

llpses  were  all  approximately  1.6  km  . The  relocation  procedure  was 
necessary  because  the  resulting  phase  velocities  are  strongly  depen- 
dent on  the  relative  event  separation.  The  PDE  and  relocated  epi- 
centers are  plotted  in  Figure  3. 

We  read  the  WWSSN  film  chips  for  all  five  events  in  our  attempt 
to  determine  whether  the  source  mechanisms  were  isotropic  and  explo- 
sive. We  found  between  21  and  29  readable  LP  polarities  for  each 
event  all  of  which  were  compressions! ; i.e.,  up.  There  was  no  evi- 
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dence  of  anomalous  amplitude  patterns  with  azimuth  indicative  of  sig- 
nificant strain  energy  release  or  of  propagating  ruptures.  Further- 
more, an  examination  of  the  long  period  horizontal  WWSSN  records  from 
a set  of  stations  approximately  40°  from  Novaya  Zemlya  and  well  dis- 
tributed in  azimuth  showed  that  no  strong  SH  waves  were  excited  by 
these  events  as  is  commonly  the  case  from  explosions  at  NTS. 

The  phase  velocity  data  from  the  six  possible  north-south  event 
pairs  are  shown  plotted  in  Figure  5.  The  error  limits  are  due  to  the 
observed  spread  in  phase  velocity  at  each  frequency.  It  appears  that 
some  of  this  spread  is  the  result  of  residual  epicenter  errors  because 
data  from  each  event  pair  tended  to  parallel  that  from  the  others. 
The  rms  error  of  the  data  over  the  band  shown  in  Figure  5 was  0.082 
km/sec.  This  figure  includes  the  bad  point  at  64  sec  period  which  our 
model  was  unable  to  fit.  Typically,  the  angle  between  the  event  and 
the  great  circle  path  to  ALPA  was  approximately  7°.  This  leads  to  a 
systematic  error  in  phase  velocity  of  -0.75Z. 

Results  of  the  Data  Inversion  Experiment 

We  inverted  the  dispersion  data  along  with  its  estimated  errors 
to  obtain  the  theoretical  dispersion  curve  indicated  by  the  solid  line 
in  Figure  5.  This  was  done  with  a generalized  inverse  program  operat- 
ing on  the  phase  velocity  partial  derivatives  (e.g.,  Wiggins,  1972; 
Rodi,  et  al.,  1975).  Our  best  fitting  model,  whose  layer  parameters 
are  given  in  Table  2,  was  obtained  by  varying  the  P and  S velocities 
of  a Gutenberg  continent  (Dorman,  et  al.,  1960)  starting  model. 
Several  tries  were  made  varying  the  crustal  thickness  in  each  case. 
The  best  fitting  model  in  terms  of  rms  error  of  fit  and  absence  of  os- 
cillations between  adjacent  layer  parameters  had  a 45  km  thick  crust. 
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Its  rms  error  of  fic  Co  che  mean  data  values  was  0.045  km/sec.  Howev- 
er, our  data  can  also  be  fit  by  a model  with  a 40  km  chick  cruse  and  a 
3X  reduedon  In  Che  shear  wave  veloclcy  of  Che  firsc  layer  below  che 
Moho.  The  rms  error  of  chls  flC  Is  slighely  greaCer  buC  still  wichln 
Che  rms  error  of  che  daca. 

The  besc  fleeing  model  Is  shown  ploeced  (on  a logarlChmic  scale) 
In  Figure  6.  A comparison  beeween  Che  S veloclcy  profiles  of  our 
model  and  Che  Gucenberg  continent  sCarcing  model  Is  shown  in  Figure  7. 
Averaging,  or  resoluCion,  kernels  for  Che  firsc  four  layer  S veloci- 
Cles  are  shown  in  Figure  8.  These  are  measures  of  che  independence  of 
Che  corresponding  parameCer.  The  firsc  Chree  S veloclCles  are  well- 
decermlned  by  chls  data  set;  however,  the  fourth  anc  , as  it  turns  out, 
succeeding  S velocities  are  relatively  poorly  resolved.  Our  best  fit- 
ting model,  therefore,  is  principally  a crust  and  uppermost  mantle 
structure.  Nevertheless,  the  data  do  not  require  an  S wave  low  velo- 
city zone  in  the  upper  mantle  in  order  to  provide  a satisfactory  fit 
for  our  preferred  model. 

Discussion 

There  is  little  information  available  in  the  literature  on  the 
crustal  structure  of  Novaya  Zemlyn.  Geologic  maps  (e«g.,  Yanshin, 
1966)  indicate  that  the  island  is  a northern  extension  of  the  Ural 
mountain  belt.  However,  whereas  the  Urals  are  a generally  linear 
feature  in  the  continental  USSR,  they  are  twisted  in  an  S-shaped  curve 
towards  Che  edge  of  the  continental  shelf.  On  Novaya  Zemlya  rocks  of 
the  western  Paleozoic  raiogeosyncline  are  exposed  (Hamilton,  1970). 
Thus,  we  would  infer  the  crustal  structure  of  Novaya  Zemlya  to  be  a 
modified  version  of  the  Ural  mountains  structure. 
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Contour  maps  of  crustal  thickness  for  the  USSR,  reproduced  as 
Figures  9 and  10,  (Kosminskaya,  et  al.,  1969;  Volvovskii,  1973)  clear- 
ly show  the  relation  of  Novaya  Zemlya  to  the  Ural  mountain  belt.  How- 
ever both  Indicate,  from  the  plotted  contours,  that  the  crust  in  the 
southernmost  island  is  between  35  and  AO  ka  thick.  The  map  by  Kosmin- 
skaya  et  al.  Includes  all  of  Novaya  Zealya  between  the  35  and  AO  ka 
thick  contours.  On  the  other  hand,  the  map  by  Volvovskii,  presumably 
the  result  of  more  analysis,  places  the  35  km  contour  near  the  fiord 
separating  the  southern  island  from  the  rest  of  Novaya  Zealya.  Howev- 
er, the  deep  seismic  sounding  profiles  indicated  by  Kosainskaya  et 
al.,  which  were  used  to  produce  these  maps,  are  all  well  to  the  south 
of  Novaya  Zemlya.  A crustal  thickness  of  50  km  is  indicated  for  the 
southern  Urals,  based  on  the  work  of  Khalevin,  et  al.  (1966).  They 
also  give  the  P velocities  just  above  and  just  below  the  Moho  as  6.4 
and  8.2  km/sec  respectively. 

The  Novaya  Zemlya  structural  model  we  present  here  (Table  2) 
agrees  in  some  respects  with  that  derived  by  Khalevin,  et  al.  (1966) 
for  the  southern  Urals.  In  particular,  our  crustal  thickness  is  A5  km 
compared  to  50  km,  our  lower  crust  P velocity  is  6.55  km/sec  compared 
to  6. A km/sec,  and  our  Moho  P velocity  is  8.18  compared  to  8.2  km/sec. 
Also,  our  third  layer  is  26  km  thick  while  theirs  is  32  km  thick.  Un- 
fortunately there  appear  to  be  no  pertinent  S velocities  available  for 
comparison  and,  since  our  model  is  primarily  an  S wave  structure, 
these  would  be  desirable  to  have. 

Conclusions 

We  have  demonstrated  a method  for  isolating  the  Rayleigh  wave 
phase  velocity  dispersion  of  a source  region  from  teleseismic  Rayleigh 
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wave  recordings  and  applied  le  to  the  southern  island  in  Novaya  Zem- 
lya.  The  preferred  crust  and  upper  mantle  structure  we  derived  by  in- 
verting our  phase  velocity  data  is  characterized  by  a 45  km  thick 
crust  and  the  absence  of  a shear  wave  low  velocity  zone  in  the  upper 
mantle*  Although  little  information  is  available  in  the  literature 
concerning  the  crustal  structure  of  Novaya  Zemlya,  our  structure  is 
compatible  published  P-wave  structures  for  the  southern  Ural  moun- 
tains. We  conclude  that  this  method  is  capable  of  providing  reliable 
Information  about  inaccessable  source  region  structures. 
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Figure  Captious 


1.  Ceometry  for  the  two  event,  single  station  method.  D is  the 
event  separation.  Deff  is  the  effective  separation  which  is  the 
projection  of  the  event  separation  on  the  great  circle  path  to 
the  station. 

2.  Map  showing  the  locations  of  the  two  test  sites  on  Novaya  Zemlya 
and  their  common  propagation  path  to  the  Alaskan  Long  Period  Ar- 
ray (ALFA). 

3.  Map  of  the  southern  part  of  Novaya  Zemlya  showing  the  PDE  and  re- 
located epicenters  of  the  five  events  used  in  this  investigation. 

4.  ALPA  center  sensor  seismograms  from  the  five  events  used  in  this 
investigation.  All  data  are  plotted  to  the  same  amplitude  and 

time  scales. 

5.  Rayleigh  wave  phase  velocity  dispersion  data,  estimated  errors, 
and  the  computed  dispersion  from  our  Novaya  Zemlya  structural 

model. 

6.  Logarithmic  scale  plot  of  the  layer  parameters  from  our  Novaya 
Zemlya  structural  model.  First  layer  parameters  are  outside  the 
range  of  the  plot. 

7.  Logarithmic  scale  plot  of  the  shear  velocity  from  our  Novaya  Zem- 
lya structural  model  compar.  to  that  from  the  starting  Gutenberg 
continent  model.  First  layt.  parameters  are  outside  the  range  of 
the  plot. 

8.  Shear  velocity  averaging  kernels  for  the  first  four  layer  shear 
velocities  as  functions  of  the  shear  and  conpressional  wave  velo- 
cities in  the  model. 

9.  Moho  depth  contour  map  of  the  USSR.  Taken  from  Kosmlnskaya,  et 
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al.,  1969. 

10.  Moho  depth  contour  nap  of  the  USSR.  Taken  from  VoLvovskli,  1973. 
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TABLE  1 
PDE  Epicenters 


Event 

Origin  Time 

Latitude 

Longitude 

“b 

8/28/72 

05  59  56.5 

73.31* 

55.09 

6.3 

11/2/7** 

OU  59  56.7 

70.82 

5l*.  06 

6.7 

8/23/75 

08  59  57-9 

73-37 

5L.6L 

6.1* 

10/18/75 

08  59  56.3 

70.8U 

53.69 

6.7 

10/21/75 

11  59  57.3 

73-35 

55.08 

6.5 

TABLE  2 


Novaya  Zemlya  Structural  Model  Layer  Parameters 


Layer 

Thickness 

(km) 

P Velocity 
(km/sec ) 

S Velocity 
(km/sec) 

Density 

(km/sec 

1 

9 

6.01 

2.83 

2.74 

2 

10 

6.13 

3.16 

2.74 

3 

2 6 

6.55 

3-98 

3.00 

4 

32 

8.18 

4.36 

3-33 

5 

30 

8.07 

4.45 

3-37 

6 

75 

7-90 

4.44 

3.42 

7 

75 

8.10 

4.45 

3.48 

8 

100 

8.50 

4.60 

3.55 

9 

100 

9-00 

4.95 

3.65 

10 

CO 

9-75 

5-37 

3-95 

DEPTH  (km) 


DEPTH  (km) 


Figure  7- 
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Summary.  A two-dimensional  finite  element  model  was  developed  for  the 
source  region  of  the  San  Fernando  earthquake.  Stochastic  inversion  of  the 
surface  displacement  data  of  Alewine  was  carried  out  to  obtain  estimates  of 
the  displacements  and  stress  drops  along  the  actual  fault  surface  in  the  finite 
element  model.  We  calculate  an  average  slip  of  222  cm  with  a rms  fit  to  the 
data  of  8 cm.  The  average  computed  stress  drop  was  290  bar,  with  a maxi- 
mum of  650  bar.  Using  these  calculated  stresses  in  a dynamic  model  of  the 
earthquake,  we  compute  theoretical  accelerograms  for  the  Pacoima  Dam  site. 
For  frequencies  less  than  2 Hz,  we  found  that  the  observed  accelerograms 
were  fitted  best  by  a model  with  a propagating  source  having  a rupture  velo- 
city of  approximately  2.5  km  s'1.  These  results  suggest  that  the  dynamic 
finite  element  method  can  be  used  to  estimate  strong  earthquake  ground 
motion  from  extended  sources  (earthquakes)  in  many  different  complex  geo- 
logic structures. 

Introduction 

The  San  Fernando,  California  earthquake  of  1971  February  9 provided  seismologists  with  a 
large  quantity  of  unusually  accurate  near-field  observational  data.  Despite  the  numerous 
investigations  to  date,  there  remains  considerable  controversy  about  its  fault  mechanism 
particularly  the  rupture  history  and  the  identification  of  individual  seismic  phases  in  the 
accelerograms.  The  present  study  is  not  intended  to  settle  the  controversy.  Rather,  we 
demonstrate  how  a two-step  procedure:  (i)  an  inversion  calculation  based  on  a static  two- 
dimensional  finite  element  method  (FEM)  with  a simple  mesh  configuration  and  faulting 
model  and,  (ii)  a dynamic  finite  element  (DFEM)  calculation  with  the  same  model,  can  be 
used  to  predict  ground  motion  for  extended  earthquake  sources  such  as  the  San  Fernando 
event.  In  constructing  best-fitting  models,  both  static  and  dynamic,  for  the  San  Fernando 
earthquake,  considerable  insight  is  gained  as  to  the  nature  of  the  faulting  process  and  the 
interpretation  of  the  energy  arrivals  recorded  by  the  Pacoima  Dam  accelerograms. 

* Formerly  at:  Geophysics  Section,  Department  of  Geosciences.  Pennsylvania  State  University,  University 
Park,  Pennsylvania  16802,  USA. 
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Choice  of  FEM  mesh  and  fault  model 

Studies  of  the  fault-plane  solutions  for  the  main  San  Fernando  event  and  the  aftershock 
sequence  (Allen  et  ai  1971 ; Wesson,  Lee  and  Cibbs  1971 ; Wyss  & Hanks  1972;  Whitcomb 
et  al.  1973;  Hanks  1974)  indicate  that  faulting  took  place  on  a three-dimensional  curved 
surface  dipping  20-25°  near  the  Earth’s  surface  and  increasing  to  50-55°  with  depth.  The 
hypocentre  is  placed  at  depths  of  8-13  km  by  these  authors.  The  maximum  width  of  the 
initial  fault,  estimated  from  the  spatial  distribution  of  the  aftershock  sequence,  is  15  km. 
The  fault  plane  solutions  such  as  Whitcomb  (1971)  obtained  (strike  N64°  W,  dip  52°  NE, 
rake  64°)  and  the  observed  ground  displacements  both  indicate  that  the  fault  mechanism 
was  predominantly  thrusting  in  a south-west  direction. 

In  order  to  make  our  finite  element  experiment  computationally  feasible,  it  is  necessary 
to  develop  a two-dimensional  model  of  this  fault  and  its  surrounding  geologic  structure. 
Jungels  & Frazier  (1973)  have  shown  that  a test  of  the  applicability  of  two-dimensional 
(plane  strain)  analysis  is  the  ratio  of  the  fault  length  to  fault  width  of  the  actual  earthquake 
being  modelled.  For  models  where  this  ratio  is  greater  than  unity,  they  found  that  plane 
strain  analysis  is  accurate  to  within  a few  per  cent,  if  we  construct  a model  with  a focal 
depth  for  the  initial  point  of  rupture  at  9.5  km,  with  a fault  plane  dipping  at  45°,  we  get 
13.4  km  for  the  fault  length  and  the  aspect  ratio  is  approximately  one.  This,  of  course,  is 
only  a geometric  criterion  as  to  the  appropriateness  of  the  plane  strain  approximation. 
However,  the  fact  that  the  fault  mechanism  indicates  a predominance  of  thrust  faulting 
suggests  that  by  considering  a profile  perpendicular  to  the  strike  of  the  faulting,  we  may 
ignore  edge  effects. 

To  construct  our  finite  element  model,  we  consider  the  geologic  structure  along  the 
profile  A-B  in  Fig.  1.  In  the  region  of  the  San  Fernando  fault,  the  profile  which  passes 
through  the  epicentre  and  the  Pacoima  Dam,  corresponds  to  the  BB'  profile  of  Jungels  & 
Frazier  (1973)  as  well  as  the  profile  used  by  Alewine  & Jordan  (1973). 

Fig.  2 shows  the  actual  finite  element  mesh  superimposed  on  the  geologic  profiles.  The 
elastic  properties  of  the  six  rock  types  were  obtained  from  well-drilling  data  published  by 
Duke  et  aL  (1971).  The  large  elements  are  placed  around  the  periphery  of  the  mesh  to  elimi- 
nate reflections  in  the  dynamic  calculations  given  below.  No  attempt  was  made,  therefore, 
to  include  mantle  structure  in  the  mesh. 

The  fault  model,  which  is  similar  to  that  used  by  Espinoza.  Harding  & Lopez-Arroyo 
(1973),  allows  the  fault  surfaces  to  move  along  planes  at  45°  to  the  Earth’s  surface.  Further- 
more, only  the  San  Fernando  fault  is  allowed  to  move;  i.e.  the  San  Gabriel  fault  remains 
locked.  Fig.  3 is  a schematic  diagram  of  the  kinematics  of  this  fault  model.  In  the  figure,  the 
direction  of  motion  of  the  two  fault  surfaces  is  indicated  by  the  unit  vectors  rq  and  nt. 

Obviously  this  is  not  a very  complicated  fault  model.  It  consists  of  only  a single  unbent 
fault  whose  plane  of  motion  is  ‘locked’  into  the  mesh,  i.e.  not  permitted  to  warp  or  translate 
relative  to  the  coordinate  system.  The  stresses  necessary  to  deform  the  mesh  into  its  final 
configuration  are  taken  to  be  the  negative  of  the  tectonic  stresses  which  were  released  by  the 
earthquake.  For  the  model,  the  stresses  are  applied  to  produce  the  deformation.  But,  for  the 
earthquake,  the  excess  tectonic  stresses  are  balanced  by  equal  and  opposite  frictional  stresses 
along  the  fault  until  they  are  removed  by  the  physical  fault  mechanism  which  produced  the 
earthquake.  Thus,  although  stress  in  our  model  is  applied  rather  than  released,  the  values  are 
equivalent.  Historically,  this  description  of  the  faulting  process  is  due  to  Reid  (republished 
1969)  and  is  known  as  the  elastic  rebound  theory.  With  the  advent  of  the  ‘new  global  tec- 
tonics’, it  is  the  generally  accepted  (simplified)  description  of  the  mechanism  of  earth- 
quakes (Brune  1974). 
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Figure  I.  ('.eographic  location  of  FEM  profile,  indicated  by  dashed  line.  Background  map  from  tirantr 


Figure  2.  Finite  element  mesh  superimposed  on  geologic  profile  shown  in  Fig.  1.  Location  of  the  Pacoima  Dam  is  indicated  by  the  letters  PD,  and  the  extent  of  the 
profile  shown  in  the  perspective  plots  is  indicated  by  dashes. 
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Figure  3.  Schematic  diagram  of  the  fault  model  showing  double  row  of  nodes  which  are  constrained  to 
move  parallel  to  the  fault . 


Static  finite  element  analysis 

Taken  together,  the  mesh  and  fault  model  provide  the  basis  of  a straightforward  static  FEM 
problem  in  the  form  of 

A'x  * f (1) 

where  A represents  the  stiffness  matrix  for  the  mesh,  which  is  itself  composed  of  individual 
element  stiffnesses  computed  from  the  static  properties  of  each  element  material  (e.g.  see 
Zienkiewicz  1971;  McCowan  1975  for  computational  details!,  x is  the  vector  of  node  point 
displacements,  and  f is  the  node  point  force  vector  (the  unknown  in  our  problem)  which 
generates  the  displacements.  Formally  we  may  rewrite  (1)  as 

/r'(=x  (2) 

and  ask:  given  a set  of  known  displacements  x,  what  are  the  corresponding  forces  f?  For 
this  static  problem,  the  appropriate  rhs  vector  is  a profile  of  surface  displacements.  Applying 
this  approach  to  the  San  Fernando  earthquake,  we  used  the  vertical  surlace  displacements 
reported  by  Alewine  (1974)  which  were  obtained  from  levelling  and  gravity  surveys. 

To  set  up  a canonical  inverse  problem,  equation  (1)  must  first  be  manipulated  in  order  to 
separate  the  known  from  unknown  quantities.  This  can  be  done  through  partitioning.  By 
identifying  the  vertical  surface  and  fault  displacements  in  the  node  displacement  vector  x 
and  the  fault  forces  in  the  force  vector  f,  the  system  can  be  written: 
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where  the  displacements  of  the  remaining  nodes  in  the  mesh  are  in  xR  and  the  forces  on  all 
the  nodes  in  the  mesh  not  on  the  fault  surface  are  zero.  This  problem  can  be  formulated  as 
an  inverse  problem  by  solving'  the  equations  to  give 


K" 


— 

B 

x 


(*») 


where  the  inverse  stiffness  matrix  (sometimes  called  the  compliance  matrix)  has  been  parti- 
tioned in  the  same  way  as  the  vectors.  The  vertical  surface  displacements  of  the  problem  are 
therefore  related  to  the  generating  force  system  acting  on  the  fault  nodes  by 

Bf,  - x,.  (5) 

The  B matrix  is  the  FEM  Green’s  function  for  the  earthquake,  i.e.  it  is  the  set  of  vertical 
displacement  profiles  for  a set  of  unit  forces  on  the  fault  nodes. 

Equation  (S)  now  has  the  three  key  elements  of  a linear  inverse  experiment.  It  has  a 
‘theory*,  expressed  in  terms  of  the  matrix  B,  which  is  to  be  used  in  explaining  the  data.  It 
has  a ‘model’  in  terms  of  the  fault  force  vector  ff  which  constitutes  the  set  ol  unknowns  in 
the  experiment,  i.e.  the  set  of  actual  values  of  the  parameters  needed  to  complete  the 
theory.  Thirdly,  it  produces  a trial  ‘data’  vector  x,  for  each  model  by  the  indicated  matrix 
multiplication  which  can  be  compared  with  the  actual  data  vector  to  judge  the  accuracy  of 
the  ‘fit’.  Expressed  in  this  form,  the  experiment  can  be  inverted  by  standard  optimization 
methods  to  yield  a description  of  the  fault  behaviour.  Because  the  problem  is  overdeter- 
mined (47  data  points  are  to  be  explained  by  25  unknowns),  the  weighted  least-squares 
(WLS)  procedure  is  the  obvious  first  approach.  It  has  the  advantage  that  additional  informa- 
tion is  introduced  into  the  estimate  by  including  the  relative  errors  of  observation  in  the 
data.  A summary  of  this  method  is  presented  in  the  Appendix. 

Fig.  4 shows  the  vertical  displacements  of  the  free  surface  calculated  along  the  profile. 
Also  shown  in  Fig.  4 are  the  observed  data  points  which  are  represented  by  circles.  The  lines 
in  the  circles  are  the  estimated  observational  errors.  This  solution  fits  the  observed  data  with 
an  rms  error  of  2.5  cm.  Fig.  5 shows  the  corresponding  slip  along  the  fault  as  a function  of 
depth.  It  is  obvious  that,  although  the  rms  error  of  the  overall  fit  is  small,  the  calculated  slip 
on  the  fault,  with  maximum  values  of  210  m,  is  quite  unrealistic.  These  slips  correspond  to 
calculated  stress  drops  of  the  order  of  tens  of  kilobars  which  are  not  plausible  in  a shallow 
earthquake. 
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Figun  4.  Static  vertical  displacement  of  the  free  surface.  Displacements  calculated  from  the  FEM  model 
using  WLS  approach  are  indicated  by  the  solid  line,  the  circles  represent  the  observed  data  from  Alewine 


(1974),  the  vertical  bars  indicate  the  uncertainty  in  the  latter. 
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Figure  S.  Relative  displacement  along  the  fault  surface  calculated  from  the  FEM  model  using  WLS 
approach. 
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The  difficulty  with  the  WLS  approach  is  that  the  theory  is  not  capable  of  accurately 
representing  the  behaviour  of  the  surface  and  fault  plane  at  the  same  time.  In  other  words, 
the  equations  are  underconstrained.  One  way  out  of  this  difficulty  (aside  from  constructing 
a whole  new  FEM  model  with  different  assumptions)  is  to  relax  the  'theory'  by  adding  a 
tradeoff  correction  term  to  the  least-squares  procedure.  The  method,  weighted  least-squares 
with  tradeoff  correction  (WLSTC),  is  a result  of  treating  the  inversion  procedure  as  a 
stochastic  process  which  allows  errors  in  the  theory  as  well  as  errors  in  the  data  when  formu- 
lating the  least -squares  equations.  In  this  form,  it  is  due  to  Foster  ( 1961)  but  it  is  also  related 
to  other  inversion  procedures  (e.g.  Backus  & Gilbert  1970;Jordan  & Franklin  1971).  Adis- 
cussion  of  this  method  is  also  included  in  the  Appendix. 

To  implement  the  stochastic  inverse  procedure,  it  is  necessary  to  determine  an  appropri- 
ate value  for  the  tradeoff  parameter  a.  From  equation  (A7),  it  is  apparent  that  a depends  on 
the  coefficient  matrix  B.  Unfortunately  no  general  way  to  set  a seems  to  work,  but  a very 
good  trial  and  error  procedure  can  be  applied.  If  the  magnitude  of  the  estimate  vector  | f | is 
plotted  against  the  rms  enor  of  the  fit 

where  M is  the  number  of  surface  displacement  observations  x,,  a curve  results  which 
resembles  a hyperbola.  This  curve,  for  the  San  Fernando  earthquake  problem,  is  shown  in 
Fig.  6.  As  Gilbert  (1970)  noted:  ‘the  place  to  be  is  down  at  the  comer’  (of  the  tradeoff 
curve).  Thus  the  ‘optimum’  value  of  a in  this  case,  a ®1  x 10"16,  produces  an  estimate  of  f 
which,  in  terms  of  slope,  is  halfway  between  the  two  extremes  of  the  curve. 


Figure  6.  Trade-off  parameter  curve.  The  value  of  the  tradeoff  parameter  a*lx  10'“  is  optimum. 
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The  results  of  using  the  WLSTC  procedure  for  this  optimum  viJue  of  the  trade-off  para- 
meter a are  shown  in  Figs  7-10.  Fig.  7 shows  the  calculated  and  observed  vertical  surface 
displacements  The  rms  error  in  fitting  the  observed  data  is  now  8.0  cm.  The  most  visible 
effect  of  using  the  optimum  tradeoff  correction  is  not  forcing  the  model  to  fit  the  two  data 
points  approximately  9 km  north-east  of  the  surface  faulting.  The  corresponding  horizontal 

O 


Figure  7.  Static  vertical  displacements  of  free  surface  calculated  using  optimum  value  of  o in  the  WLSTC 
procedure. 
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Figure  8.  Horizontal  displacements  of  free  surface  corresponding  to  those  in  Kig  7. 


Figure  10.  Calculated  itren  drop  corresponding  to  displacements  shown  in  Fig.  9.  The  dashed  line  indi- 
cates anomalous  region  of  FEM  stress  calculation. 
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surface  displacements  are  shown  in  Fig  8.  Mere  the  use  of  the  tradeoff  parameter  results 
in  a displacement  profile  closely  resembling  that  of  an  ideal  thrust  fault  The  calculated  dis- 
placements along  the  fault  surface  are  shown  in  Fig.  9.  The  maximum  slip  is  now  4.5  m at  a 
point  1.56  km  down  the  fault,  with  a second  local  maximum  of  3.7  m located  at  7.23  km 
down  the  fault  The  maximum  computed  stress  drop  (Fig  10)  is  b50  bar,  at  a point  2.4  km 
from  the  free  surface.  The  average  stress  drop  is  290  bar  (The  apparent  peak  in  the  stress 
drop  at  the  surface,  indicated  by  the  dashed  part  of  the  curve  in  Fig  10,  should  be  disre- 
garded as  it  is  an  artefact  of  the  way  stress  is  computed  from  (he  spatial  derivatives  of  (he 
displacement  solution.) 

As  a check,  we  recomputed  an  inverse  solution  using  a tradeoff  parameter  equal  to  one- 
fifth  of  the  optimum  value  The  mis  lit  of  this  solution  was  5.7  cm.  The  fault  slip  and  stress 
drop  which  resulted  from  the  non-optimum  tradeoff-corrected  model  follow  the  overall 
behaviour  of  the  corresponding  optimum  solutions  Thus,  over  the  tange  shown  in  these 
examples,  the  principal  results  of  the  W'LSTC  models  are  relatively  insensitive  to  a change 
in  the  tradeoff  parameter  of  a factor  of  5 in  the  vicinity  of  the  'optimum’  value. 


Dynamic  finite  element  analysis 

If  the  mesh  does  not  change  with  time,  then  the  dynamical  problem  mesh  can  be  developed 
from  the  static  case  by  simply  assigning  mass  and  damping  to  each  element  These  quantities, 
however,  must  be  distributed  throughout  each  element  in  a manner  consistent  with  the  stiff- 
ness distribution  (McCowan  1975)  The  resulting  set  of  equations  is 

Ma+Dv  + A x - f 

where  M and  D are  the  mass  and  damping  matrices  and  a and  v are  the  node  point  accelera- 
tion and  velocity  vectors,  respectively  To  complete  the  definition  of  the  problem,  relations 
linking  acceleration  and  velocity  to  displacement  must  be  supplied.  A common  set  of  such 
relations  are  the  finite  difference  equations  due  to  Newmark  (1950) 

h 

»«.i  " +-  +*n.t) 

X*-,  » x„  +\nh  + (j  -d!A:an  *3*Ja„  * , 

where  *1  is  the  integration  parameter,  and  h the  integration  time  step  With  these  equations 
a,  v.  and  x at  time  n ♦ 1 can  be  computed  from  their  values  at  time  n.  Throughout,  a value 
^■0.25  was  used  as  this  leads  to  an  unconditionally  stable  integration  scheme  (McCowan 
1975).  The  integration  time  step  h was  chosen  so  as  not  to  exceed  the  Courant  limit  given  by 

J 

h < - 
a 

where  a is  P- wave  velocity,  and  J is  the  smallest  mesh  spacing  A damping  matrix  appiopnate 
to  this  integration  scheme  and  suitable  for  use  in  problems  concerned  with  Rayleigh  waves 
(Frazier  et  aL  1973;  McCowan  1975)  is 

h 

D - - A 

2 

In  this  form,  the  damping  matrix  attenuates  waves  proportional  to  their  frequency  squared 
and  its  effects  are  most  pronounced  in  coarse  regions  of  the  mesh  Since  the  large  elements 
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in  Fig.  2 are  placed  around  the  periphery  of  the  mesh  to  prevent  reflections  from  the  bound- 
aries occurring  during  the  time  window  of  interest,  the  effect  of  damping  is  minimal  in  the 
region  of  the  Pacoima  Dam. 

The  optimum  fault  force  model  derived  from  the  static  calculation  above  provides 
important  constraints  on  the  dynamics  of  the  San  Fernando  earthquake.  Specifically,  it 
limits  the  number  of  possible  dynamic  fault  models  by  setting  the  asymptotic  behaviour  of 
the  faulting  process.  That  is,  the  force  time  history  acting  on  each  fault  node  must  approach 
the  calculated  static  value  for  large  times.  The  time  history  itself  is  therefore  an  unknown  in 
the  dynamic  inversion  experiment.  Unfortunately,  the  number  of  possible  time  histories  is 
too  large  to  be  economically  inverted  by  formal  methods  so  the  problem  was  parameterized 
by  assuming  various  functional  forms  for  the  stress  time  history  and  comparing  the  observed 
and  calculated  accelerograms. 

The  simplest  earthquake  stress  time  history  is  a propagating  step  function  of  applied 
stress.  Following  Hanks  (1974),  the  initial  rupture  is  placed  in  the  hypocentral  region  and  is 
assumed  to  propagate  along  the  fault  towards  the  surface.  As  the  rupture  reaches  each  fault 
node,  the  applied  force  goes  to  the  static  or  asymptotic  value  instantaneously.  Symbolically, 
this  source-time  function  is 

f (s,  r)  = f„  (s)  H(t  - s/VR) 

where  s is  the  distance  along  the  fault  from  the  hypocentre,  f„  is  the  static  force,  H is  the 
unit  Heaviside  step  function,  and  VR  is  the  rupture  velocity.  In  this  idealized  earthquake 
mechanism  model,  the  only  undetermined  parameter  is  the  rupture  velocity. 

The  results  of  trying  two  rupture  velocities,  2.0  and  2.5  km  s'1  are  presented  in  Figs 
11-17.  Fig.  11  is  a pair  of  perspective  plots  of  the  horizontal  and  vertical  components  of 
surface  displacement  along  the  profile  indicated  by  dashes  in  Fig.  2 for  the  2.5  km  s'1  rup- 
ture velocity.  The  corresponding  plots  for  the  2.0-km  s~‘  case  (not  shown)  are  similar.  It 
appears  that  most  of  the  energy  intersects  the  surface  as  a vertically  travelling  wave  because 
the  apparent  velocity  of  the  principal  disturbance  on  the  surface  (>  10  km  s’1)  is  higher  than 
the  P- wave  velocity  in  any  element  of  the  mesh.  There  is  no  evidence  of  significant  energy 
propagating  along  the  Earth’s  surface  from  the  fault  break  towards  the  Pacoima  Dam.  The 
perspective  plots  also  show  that  the  surface  motion  has  a local  maximum  near  Pacoima  Dam. 
Ground  motion  on  either  side  of  this  peak  is  about  a third  smaller  than  at  the  Pacoima  Dam 
node. 

Fig.  12  is  a schematic  diagram  of  possible  energy  paths  from  the  hypocentre  to  the 
Pacoima  Dam.  The  travel  times  1 1 and  f2  are  for  the  direct  P and  S waves.  Times  r3  and  r4 
include  the  time  of  the  rupture  propagating  up  the  fault  and  the  time  of  seismic  energy  to 
propagate  from  the  fault  to  the  Pacoima  Dam.  For  the  1 3 path,  the  seismic  energy  is  assumed 
to  propagate  as  a shear  wave  from  the  point  on  the  fault  nearest  the  Pacoima  Dam.  The  r4 
path  is  the  time  of  a Rayleigh  wave  to  reach  the  receiver  from  the  fault  after  the  rupture  has 
broken  the  surface.  These  times  are  indicated  on  the  figures  which  follow  showing  the 
calculated  ground  motion  at  the  Pacoima  Dam  location. 

The  principal  results  of  the  dynamic  experiment  are  shown  in  Figs  13—16.  Both  compon- 
ents of  acceleration,  velocity,  and  displacement  are  plotted  for  each  of  the  two  rupture 
velocities.  The  major  arrival  for  both  cases  occurs  around  r3.  The  direct  P and  5 body-wave 
phases  are  smaller  than  the  f3  arrival  by  an  order  of  magnitude  in  both  cases.  Furthermore, 
the  result  noticed  on  the  perspective  plots  is  confirmed:  there  are  no  large  amplitude  phases 
at  or  later  than  1 3.  The  particle  motion  of  the  r3  part  of  the  time  history  is  up  and  south  for 
both  rupture  velocities,  i.e.  parallel  to  that  of  the  hanging  wall  side  of  the  fault.  This  is  con- 
sistent with  the  arrival  being  due  to  a shear  wave  propagating  up  and  north  in  a direction 
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Kigur*  1 1.  Calculated  dynamic  ground  displacements  at  the  free  surface  plotted  at  a function  of  distance 
north  from  the  fault,  for  a rupture  velocity  of  2.5  kms'V  PD  indicates  location  of  Pi  count  Dam  Positive 
V and  IV  amplitude!  indicate  displacement  northward  and  upwards,  respectively 
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Figure  12.  Possible  energy  paths  from  hypo  cent  re  to  Pacoima  Dam. 


perpendicular  (o  (he  fault.  Clearly  then,  this  source  model  favours  the  t}  arrival  hypothesis: 
the  principal  energy  reaching  Pacoima  Dam  is  the  result  of  the  shear  waves  generated  by  the 
rupturing  in  the  region  of  the  fault  nearest  the  receiver. 

A comparison  of  the  computed  accelerograms  with  those  observed  at  Pacoima  Dam  is 
presented  in  Fig.  17.  The  observed  accelerograms  have  been  low-pass  filtered  with  a phase- 
less  filter  having  a corner  frequency  of  2 Hz  and  a rolloff  of  1 2 dB/octave.  This  was  done  in 
order  to  limit  the  frequency  content  of  the  observed  data  to  that  of  the  DFEM  calculations, 
the  latter  being  restricted  by  the  resolution  of  the  mesh.  In  addition,  the  low-pass  filtering 
removes  those  frequencies  from  the  observed  accelerograms  most  likely  to  have  been 
affected  by  topography  at  the  receiver  site  (Boore  1973;  Wong  & Jennings  1975). 

The  Pacoima  Dam  accelerograms  seem  to  have  been  triggered  by  the  direct  P arrival,  a 
result  indicated  by  the  existence  of  high  frequency  energy  early  in  the  record  (Hanks  1974). 
The  computed  accelerograms,  however,  begin  at  the  initial  rupture  time.  So.  if  Hanks’  infer- 
ence is  correct,  the  traces  can  be  properly  aligned  by  subtracting  the  P travel  time  (about  2 s) 
from  our  computed  accelerograms.  Doing  this  brings  the  large  arrival  occurring  at  5 s on  (he 
Pacoima  Dam  accelerograms  into  correspondence  with  the  phase  for  the  2.5  km  s'1  rup- 
ture velocity.  The  peak  to  peak  amplitudes  agree  with  18  and  40  per  cent  for  the  vertical 
and  horizontal  respectively.  The  2.0  km  s'1  rupture  velocity  case  is  characterized  by  a r3 
phase  arriving  a second  or  two  later  and  agreeing  in  amplitude  within  33  and  8 per  cent  for 
the  vertical  and  horizontal.  These  amplitudes  are  generally  within  the  accuracy  expected  of 
the  method  (McCowan  1975).  The  timing  of  the  t3  arrival,  though,  tends  to  support  the 
higher  rupture  velocity  hypothesis.  Based  on  these  results  we  can  infer  that  the  stress  step 
assumption  is  acceptable  for  representing  the  dynamic  response  at  frequencies  below  2 Hz 
and  that  the  method  can  be  used  to  predict  approximately  the  spatial  and  temporal  distri- 
bution of  ground  motion  in  complicated  geologic  settings. 

Discussion 

It  is  instructive  to  compare  our  static  FEM  results  with  those  obtained  by  previous  investi- 
gators. Several  authors  examined  the  San  Fernando  earthquake  problems  by  inverting 
analytic  or  piecewise-analytic  solutions.  All  employed  three-dimensional  fault  models  in 
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FiguM  14.  Computed  vertical  components  corresponding  to  those  of  Fig.  13.  Arrow  indicates  motion 
into  the  ground. 
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homogeneous  media.  Canitez  & Toksoz  (1972)  and  Alewine  & Jordan  (1973)  found  static 
solutions  in  a halfspace  which  were  derived  from  surface  displacement  data.  Both  models 
produced  average  fault  slips  of  about  3 m.  The  rms  error  of  Alewine  & Jordan’s  fit  was  6 cm. 
Boore  A Zoback  (1974)  and  Trifunac  (1974)  derived  dynamic  models  from  the  strong- 
motion  accelerograph  data.  Boore  & Zoback  fitted  the  peak  velocity  in  the  data  with  a dyna- 
mic model  which  agreed  with  that  of  Alewine  & Jordan  in  the  static  limit.  Trifunac  used  a 
wholespace  solution  which  he  subsequently  doubled  to  simulate  the  effect  of  the  free 
surface.  He  obtained  a higher  average  slip  of  8-10  m.  In  contrast,  Jungels  & Frazier  (1973) 
inverted  the  observed  static  surface  displacement  data  with  a FEM  model  by  direct  methods, 
i.e.  ‘trial  and  error’.  Their  model  gave  an  average  fault  slip  much  the  same  as  that  obtained 
by  Alewine  & Jordan.  Thus,  with  the  exception  of  the  work  by  Trifunac,  our  results  for  the 
relative  displacements  along  the  fault  plane,  2.2  m on  average,  are  in  excellent  agreement 
with  those  obtained  by  previous  investigators,  as  is  the  rms  error  of  our  fit;  8.0  cm  for  the 
optimum  solution. 

The  slip  and  fault  dimensions  can  be  used  to  compute  the  seismic  moment 

Mt,  = uSA 

where  u is  the  shear  modulus  in  the  fault  zone,  S is  the  average  slip  or  displacement  across 
the  fault,  and  A is  the  faulted  area.  Since  our  model  is  two-dimensional,  and  because  both  S 
and  n vary  along  the  fault,  we  may  compute  the  seismic  moment  per  kilometre  of  fault 
length 

/•swfice 

M = 10s  I n(x)S(x)dx 

J hypo  centre 

where  /j  is  in  units  of  dyne  cm'2,  S and  x are  in  cm.  From  the  slips  given  in  Fig.  9 and  the 
seismic  parameters  of  our  model  (Fig.  2),  we  obtain  M = 0.27  x 1025dyne  cm  per  km  of  fault 
width,  which  is  a factor  of  3 less  than  the  value  obtained  by  Jungels  & Frazier  (1973)  from 
their  profile  BB'.  If  we  assume  the  faulted  area  to  be  15  km  wide,  than  the  seismic  moment 
given  by  our  results  is  Mq  = 4.0  x 1025 dyne  cm.  This  value  is  in  general  agreement  with  the 
values  obtained  by  other  investigators:  Wyss  & Hanks  (1972),  Alewine  (1974),  and  Hanks 
(1974),  for  example. 

Estimates  for  the  stress  drop  during  the  San  Fernando  earthquake  vary.  Wyss  & Hanks 
(1972)  calculated  the  stress  drop  to  be  14  bar  from  teleseismic  body-wave  data.  Canitez  & 
Toksoz  (1972)  estimated  the  stress  drop  to  be  70  bar.  Jungels  & Frazier  (1972)  estimated 
the  stress  drop  to  be  approximately  25  bar.  Alewine  & Jordan  (1974)  did  not  calculate  the 
average  stress  drop,  but  more  recent  work  (Alewine,  private  communication)  indicates  that 
it  is  on  the  order  of  100-200  bar.  In  a reinterpretation  of  the  Pacoima  Dam  accelerograms. 
Hanks  (1974)  has  postulated  massive  stress  drops  of  350-1400  bar  localized  in  the  hypo- 
central  region.  Our  static  results,  computed  directly  from  the  finite  element  model,  give  an 
average  stress  drop  of  290  bar  with  a maximum  of  650  bar  and  support  the  higher  values 
reported  by  Alewine  &.  Hanks.  However,  in  our  model,  the  highest  stress  drop  is  not  located 
at  the  initial  point  of  rupture. 

The  higher  stress  drops  are  also  supported  by  the  dynamic  results.  Because  the  FEM  prob- 
lem is  linear,  an  average  stress  drop  of  30  bar,  would  produce,  for  any  reasonable  rupture 
velocity,  maximum  computed  accelerations  one  tenth  the  value  of  those  observed.  The 
dynamic  models  predict  that  the  major,  low-frequency  energy  in  the  early  part  of  the 
Pacoima  Dam  accelerograms  is  due  to  shear  waves  generated  as  the  rupture  propagates  up 
the  fault  surface  to  points  nearest  the  receiver  at  Pacoima  Dam.  Although  the  rupture  radi- 
ates all  along  the  fault  surface,  the  peak  amplitude  in  the  accelerograms  predicted  by  this 
model  is  due  to  S-wave  energy  propagating  from  the  point  on  the  fault  nearest  the  receiver. 
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The  predicted  amplitudes  are  generally  in  agreement  with  those  observed  for  the  correspond- 
ing portion  of  the  Pacoima  records  and  the  calculated  vs  observed  arrival  times  indicate  a 
rupture  velocity  nearer  to  2.5  km  s'1  than  2.0  km  s'1. 

The  t j arrival  hypothesis  advanced  by  this  DFEM  model  agrees  with  the  results  presented 
by  Bolt  (1972)  who  found  that  different  frequency  bands  in  the  Pacoima  Dam  accelero- 
grams contained  different  energy  arrivals.  In  particular,  energy  originating  from  the  rupture 
as  it  propagated  up  the  fault  predominates  below  2 Hz.  Thus  our  dynamic  model,  con- 
strained to  low  frequencies  by  virtue  of  its  limited  bandwidth,  corresponds  to  and  agrees 
with  Bolt  s interpretation  of  the  data. 

The  local  maximum  in  the  surface  motion  near  Pacoima  Dam  shown  in  Fig.  1 1 appears 
in  similar  halfspace  problems  run  by  Geller  & Frarier  (197t>).  In  our  model,  the  effect  may 
have  been  accentuated  by  the  high  impedance  contrast  across  the  nearby  geologic  contact 
between  the  sediments  of  the  Modelo  formation  and  the  crystalline  rocks  (granite)  to  the 
north  (see  Fig.  2).  Geller  & Frazier  also  found  the  effect  of  crustal  structure  on  their  com- 
puted waveforms  to  be  significant. 

All  these  results  are,  of  course,  dependent  on  the  FEM  mesh  design,  its  resolution,  and 
number  of  degrees  of  freedom,  as  well  as  the  assumptions  implicit  in  the  dynamics  of  the 
fault  model.  In  the  region  of  the  fault  and  Pacoima  Dam.  the  average  mesh  spacing  was  on 
the  order  of  10  nodes  per  wavelength  of  the  body  waves  being  propagated.  Tins  figure  is 
commonly  used  as  a criterion  for  accuracy  in  DFEM  calculations  (e.g  NlcCowan  1975; 
Geller  & Frazier  197p).  A finer  mesh  would  immediately  serve  to  extend  the  frequency  con- 
tent of  the  arrivals  beyond  the  present  2-Hz  limit. 

The  effects  of  mesh  coarseness  on  static  FEM  calculations  are  more  difficult  to  assess. 
There  is  much  civil  engineering  literature  and  lore  on  the  subject  (e  g.  Zienkiewicz  1971)  As 
a general  rule,  coarser  elements  exhibit  more  stiffness  in  linear  elastic  problems  Pius  the 
unknown  force  vector  in  equation  (5)  and  the  corresponding  stresses  in  Fig  10  may  be  arti- 
ficially large  in  our  model.  How  much  so  can  be  estimated  by  noting,  as  above,  that  our  fault 
slip  solution  generally  agrees  in  magnitude  with  that  from  other  investigators.  Then  the 
average  strain  from  7 to  13  km  down  the  fault  is  approximately  5 x 10"*  (see  Fig  9).  Using  a 
shear  modulus  for  granite  appropriate  to  our  structure.  q = 4x  10"  dyne  cm'1,  gives  an 
average  stress  difference  of  200  bar  in  this  region  of  the  fault.  Thus  to  believe  fault  slips  as 
we  and  others  predict,  means  the  stresses  must  be  on  the  order  of  hundreds  of  bars.  Any 
improvement  due  to  tnesh  refinement  cannot  be  substantial  when  the  fault  oftset  and 
stresses  agree  as  they  do. 

Against  this  must  be  balanced  the  effects  of  our  two-dimensional  model  Since  the  profile 
shown  in  Fig.  1 passes  through  the  earthquake  epicentre,  our  computed  stresses  in  Fig  10 
will  be  the  maximum  occurring  anywhere  on  the  fault  surface  Stresses  averaged  over  ti  e 
horizontal  extent  of  the  fault  must  be  somewhat  sniallei 

The  results  of  this  FEM  experiment,  with  us  small  mesh  and  idealized  fault  model, 
demonstrate  that  it  is  feasible  to  predict  the  spatial  and  temporal  distribution  of  ground 
motion  for  earthquakes  in  other  complicated  geologic  settings  with  a reasonable  computa- 
tional effort.  While  a similar  treatment  of  three-dimensional  heterogeneous  st  me  lures  can  be 
readily  developed,  the  computational  effort  involved  is  prohibitive  except  on  the  laigest.  and 
fastest  digital  computers  presently  available 


Acknowledgments 

It  is  a pleasure  to  acknowledge  the  assistance  afforded  us  b\  R V'  Alewme  S T Maiding 
kindly  supplied  us  with  a preliminary  copy  ol  Ins  FFM  mesh  from  which  the  inesli  used  in 


184  D.  W McCowan,  P.  Clover  and  S.  S.  Alexander 

this  paper  was  developed.  Tins  research  was  supported  by  the  National  Oceanic  a. 
spheric  Administration  under  Grant  N-22-125-72(G)  and  the  US  Advanced  Research  Pro, 
Agency  under  Grant  AFOSR-73-2S15.  One  of  the  authors  (DWM)  was  also  supported  by  the 
Advanced  Research  Projects  Agency  of  the  Department  of  Defense. 


References 

Ale  wine,  R.  W.,  1974.  Application  of  linear  inversion  theory  toward  the  estimation  of  seismic  source 
parameters,  PhD  thesu.  303  pp.,  California  Institute  of  Technology,  Pasadena. 

Alewma.  R.  W.  A Jordan.  T.  H.,  1973.  Generalized  inversion  of  earthquake  static  displacement  fields, 
Geophys.  J.  R ostr  Soc..  35,  357-361. 

Allen.  C.  R.,  Engen,  G.  R.,  Hanks,  T.  C.,  Nordquist,  J.  M.  & Thatcher,  W.  R„  1971.  Main  shock  and  larger 
aftershocks  of  the  San  Fernando  earthquake,  February  9,  through  March  1,  US  geol.  Surv  Prof 
Paper  733,  17-20. 

Backus,  G.  A Gilbert,  K.,  1970.  Uniqueness  in  the  inversion  of  inaccurate  gross  Earth  data.  Phil  Trans  R 
Soc  Land.  A..  266,  123-192. 

Bolt,  B.  A.,  1972.  San  Fernando  rupture  mechanism  and  the  Pacoima  strong-motion  record.  Bull,  seism. 
Soc.  An.  62,  1053-1961. 

Boo  re,  0.  M..  1973.  The  effect  of  simple  topography:  implications  for  the  accelerations  recorded  at 
Pacoima  Dam,  San  Fernando  Valley,  California,  Buff,  seism  Soc  Am..  63,  1603-1609. 

Boo  re,  D.  M.  A Zoback,  M.  D„  1974.  Two-dimensional  kinematic  fault  modeling  of  the  Pacoima  Dam 
strong-motion  recordings  of  the  February  9,  1971,  San  Fernando  earthquake.  Bull,  seism  Soc. 
An.  64,555-570. 

Brune,  J.  N , 1974.  Current  status  of  understanding  quasi-permanent  fields  associated  with  earthquakes, 
Trans  Am  geophys.  Union,  55,  820-827. 

Canitez,  N.  A Toksoz,  M.  N.,  1972.  Static  and  dynamic  study  of  earthquake  source  mechanism.  Son 
Fernando  earthquake,/ geophys  Res.  77,2583-2594. 

Duke,  C.  M„  Johnson,  J.  A.,  Kharraz.  Y..  Campbell,  K.  W.  A Malpiede,  N.  A..  1971 . Subsurface  site  con- 
ditions and  geology  in  the  San  Fernando  earthquake  area,  UCLA  Engineering  Report.  188  pp., 
UCLA-ENG-7206,  University  of  California,  Los  Angeles. 

Espinosa,  A.  F.,  Harding,  S.  T.  A Lopez-Arroyo,  A..  1973.  Strong  motion  accelerations  and  displacements 
at  near  and  intermediate  distances;  San  Fernando  earthqua>  \ Earthquake  Notes.  XL1V,  49. 

Foster,  M..  1961.  An  application  of  the  Wiener-Kiolmogorov  sm  ihing  theory  to  matrix  inversion,  J Soc 
indust  Math..  9,387-392. 

Geller,  R.  J.  and  Frazier,  G.  A.,  1976.  Near  Held  modelling  of  dislocations  in  a heterogeneous  crust:  a 
dynamic  finite  element  approach,  / geophys.  Res.  in  press. 

Gilbert,  F.,  1970.  Inverse  problems  for  the  Earth's  normal  modes.  The  nature  of  the  solid  earth.  McGraw- 
Hill,  New  York. 

Grantz,  A.,  1971.  The  San  Fernando,  California,  earthquake  of  February  9,  1971,  US  geoL  Surv.  Prof 
Paper.  733,  1 -4. 

Hanks,  T.  C,  1974.  The  faulting  mechanism  of  the  San  Fernando  earthquake,/,  geophys  Res.  1215- 
1229. 

Jordan,  T.  H.  A Franklin.  J.  N..  1971.  Optimal  solutions  to  a linear  inverse  problem  m geophysics,  Pmc. 
Nat.  Acad  Scl  Am.  68,  291  -293. 

Jungels,  P.  H.  A Frazier,  G.  A.,  1973.  Finite  element  analysis  of  the  residual  displacements  for  an  earth- 
quake rupture:  Source  parameters  for  the  San  Fernando  earthquake,  /.  geophys  Res.  78,  5062- 
5083. 

McCowan.  D.  W.,  1975.  Dynamic  finite  element  analysis  with  applications  to  seismological  problems,  PhD 
thesis.  188  pp..  The  Pennsylvania  State  University,  University  Park.  Pasadena. 

Newmark,  N.  M.,  1959.  A method  of  computation  of  structural  dynamics,/,  eng.  mech.  Div..  Prvv.  .4m. 
Soc.  civ.  eng.,  85,67-94. 

Reid,  H.  F.,  1969.  The  California  Earthquake  of  April  IS.  190ti,  Carnegie  Institution  of  Washington. 
Washington.  DC. 

Trifunac,  M.  D.,  1974.  A three-dimensional  distortion  model  for  the  San  Fernando,  California,  earth- 
quake of  February  9,  1971, Bull  seism.  Soc.  Am.  64,149-172. 

Wesson,  R.  L.,  Lee,  W.  H.  K.  A Gibbs,  J.  F.,  1971.  Aftershocks  of  the  earthquake.  US  geoL  Surv  Pnrf. 
Paper.  733.  24-29. 


4 


Analysis  of  the  IV  71  California  earthquake  185 

Whitcomb,  J H.,  Allen,  C-  R..  Germany.  J.  D.  i Hileman.  J.  A.,  1973.  The  1971  Sen  Fernando  earth- 
quake series:  focal  mechanisms  and  tectonics,  Rev  Geophys  space  Phyt,  1 1, 693-730. 

Wong,  H.  L 3l  Jennings.  P.  C..  197}.  Effects  of  canyon  topography  on  strong  ground  motion.  Bull  seism. 
Soc  Am..  65.  1239-1257. 

Wyss,  M.  and  Hanks,  T.  C.  1972.  The  source  parameters  of  the  San  Fernando  earthquake  inferred  from 
teleseismic  body  waves, Bull  seism.  Soc  Am.  62,591-602. 

Zienkiewicz,  O.  C..  1971.  The  finite  element  method  in  engineering  science.  McGraw-Hill.  London. 

Appendix 

Least-squares  method  with  tradeoff  correction 

The  usual  least -squares  procedure  for  estimating  the  solution  of  the  M by  .\(M  > ,V)  system 
of  linear  equations 

flf  ■ x + e (Al) 

where  B is  the  M by  N coefficient  matrix,  f is  the  vector  of  estimates,  x is  the  vector  of 
observations,  and  e is  a vector  of  errors  in  each  observation,  consists  of  minimizing  the  mean 
squared  error  in  the  system 

- (eTe)  = - (xT  - fTflT ) (flf  - x).  (A2) 

At  At 

The  minimization  leads  to  the  familiar  equation  for  the  estimate  off 
f *(BTBy'BTx  (A3) 

provided  that  BJB  is  non-singular. 

A simple  modification  of  this  procedure  consists  of  weighting  the  contribution  of  each 
equation  to  the  mean  squared  error  by  its  own  error  of  observation.  Thus,  if 

E * <«eT  > (A4) 

is  the  variance-covariance  matrix,  the  mean  squared  error  becomes 

- (xT  - fTfiT)  (flf  - x).  ( A5) 

M 

The  stochastic  inverse  of  the  matrix  B defined  by  Foster  (l^ol)  can  be  thought  of  as 
that  due  to  adding  a ‘tradeoff  term  (a)  into  equation  (A5)  as  follows 

- (xT  - f rBT)  £**'  (tff  - x)  + - afT  f.  (Ab) 

M M 

Thus  a controls  how  much  of  the  additional  term  the  mean  squared  error  ‘sees'.  Minimiza- 
tion of  this  expression  gives 

f * (BTE~l  B + aiy1  BJE~l  x.  (A7) 

If  the  errors  in  the  observations  are  uncorrelated,  then  the  matrix  E is  diagonal  and  the  non- 
zero elements  are  simply  the  weights  in  the  estimate. 
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A FAST,  ACCURATE  METHOD  FOR  COMPUTING  GROUP-VELOCITY 
PARTIAL  DERIVATIVES  FOR  RAYLEIGH  AND  LOVE  MODES 

Bv  W.  L.  Rodi,  P.  Glover,  T M.  C.  Li,  and  S S.  Alexander 


ABSI  RACT 


A method  for  quickly  and  accurately  calculating  Rayleigh-  and  Love-mode 
group-velocity  partial  derivatives  with  respect  to  model  parameters  (m)  is 

„ dC 

developed.  The  method  requires  computer  codes  that  calculate  (,  t , and  — 

< m <„ 

<>C|  SU  | 

and  emplovs  numerical  differentiation  of  j to  yield  . The  method  is 

cm  |„  Cm  |„ 

CC I CU I 

fast  because  and  for  all  the  model  parameters  can  be  obtained  at  a 

< „ Cm 

given  frequency  from  only  two  solutions  of  the  period  equation.  The  accuracy 
of  the  method  is  established  with  two  examples.  For  Love  waves,  the  group- 
velocity  partials  computed  by  this  method  agree  exactly  with  those  obtained 
analytically  by  Novotny  (1970).  For  Rayleigh  waves,  comparison  with  a “brute 
force"  calculation  of  group-velocity  partials  showed  agreement  to  the  order  of 
0.0C002.  Systematic  inversion  of  group-velocity  data  separately  or  in  combination 
with  phase-velocity  data  is  computationally  feasible  using  this  method. 


Introduction 

The  development  of  linear  inverse  theory  (for  example,  Foster.  1961 ; Franklin,  1970; 
Backus  and  Gilbert,  1970;  Wiggins,  1972;  Jordan.  1973)  has  made  it  possible  to  deduce 
earth  structures  from  a variety  of  seismological  data.  The  successful  application  of  this 
theory  requires  a computationally  efficient  method  for  repeatedly  generating  partial 
derivatives  of  the  data  predicted  by  an  earth  model  with  respect  to  the  parameters  of  that 
model.  This  paper  offers  an  efficient  method  for  computing  partial  derivatives  of 
Rayleigh-  and  Love-wave  group  velocity  for  a horizontally  stratified  earth  model.  The 
method  is  accurate  and  readily  implemented  using  any  existing  computer  program  that 
calculates  phase  velocity,  group  velocity,  and  phase-velocity  partials  with  respect  to 
model  parameters. 

Compared  to  phase-velocity  partials,  group-velocity  partials  have  proved  cumbersome 
and  time-consuming  to  compute.  Also,  the  most  efficient  of  the  presently  available 
techniques  to  compute  group-velocity  partials  applies  to  Love  waves  but  not 
Rayleigh  waves.  This  has  prevented  the  calculation  of  Rayleigh-wave  group-velocity 
partials  on  a routine  basis.  A brief  review  of  the  available  methods,  presented  below  in 
order  of  increasing  efficiency,  will  point  out  these  difficulties. 

First,  the  “brute  force"  approximation  to  group-velocity  partial  derivatives  can  be 
obtained  by  calculating  the  changes  in  group  velocity  at  each  frequency  caused  by  the 
separate  perturbation  of  each  model  parameter.  This  method  can  be  used  for  both  Love 
and  Rayleigh  modes  but  is  very  time-consuming. 

Second,  Novotny  (1970)  derived  exact  expressions  for  Love-wave  group-velocity 
partials  for  a horizontally  layered  earth  model.  His  expressions  involve  first  and  second 
derivatives  of  the  root  equation  with  respect  to  phase  velocity,  frequency,  and  the  model 


U05 


w.  L.  RODI,  P.  GLOVER,  T.  M.  C.  LI,  AND  S.  S.  ALEXANDER 


parameters.  Novotny  evaluates  these  derivatives  in  terms  of  first  and  second  derivatives 
of  the  Thomson-Haskell  layer  matrices  (Haskell,  1953).  The  computer  time  required  to 
generate  all  the  necessary  layer  matrix  derivatives  and  to  perform  the  many  matrix 
multiplications  needed  is  considerable. 

Third.  Harkrider  ( 1968)  and  Anderson  and  Harkrider  ( 1968)  derived  exact  expressions 
for  Love-wave  group-velocity  partials  in  a different  way.  By  using  the  variational  principle, 
they  were  able  to  express  the  partials  in  terms  of  first  derivatives  of  energy  integrals. 
Because  only  first  derivatives  of  the  layer  matrices  are  required,  their  method  requires 
fewer  calculations  than  Novotny's.  This,  to  our  knowledge,  has  not  been  applied  to 
Rayleigh  waves. 

Finallv,  a method  due  to  Kosloff  (1975)  for  approximating  Love-wave  group-velocity 
partials  reduces  the  number  of  calculations  even  further.  Using  second-order  perturbation, 
Kosloff  (1975)  expresses  the  partials  in  terms  of  a sum  of  energy  integrals  for  all  of  the 
modes  at  a given  wave  number  or  frequency.  He  employs  a scheme  to  eliminate  leaky 
modes  and  force  the  sum  to  converge  in  five  or  fewer  terms;  thus  the  method  requires 
approximately  four  or  five  root  calculations  per  frequency  (Kosloff,  1975). 

The  method  for  obtaining  group-velocity  partials  described  in  this  paper  works  for 
any  Rayleigh-  or  Love-wave  mode.  It  requires  only  the  calculation  of  a double  set  of 
roots  and  phase-velocity  partials  to  obtain  group-velocity  partials  for  all  model  para- 
meters at  a given  frequency.  The  basis  for  the  method  is  the  theoretical  relationship 
between  group-  and  phase-velocity  partial  derivatives  which  is  derived  in  the  next 
section. 


Theory 

Define  U to  be  the  group  velocity  and  C to  be  the  phase  velocity  at  a frequency  to. 

ell  rC 

Let  and  , , respectively,  be  the  derivatives  with  respect  to  a model  parameter  m. 

c m c m „ 

holding  to  fixed,  of  group  and  phase  velocity. 

By  ignoring  all  the  model  parameters  of  a layered  half-space  except  m,  C can  be 
written  as  a function  of  to  and  m 

C=f(w,m)  (I) 

where  C is  a root  of  the  period  equation 

F(C,  co,  m)  = 0.  (2) 

If  the  function  f2  is  defined  to  be  the  partial  derivative  of/ with  respect  to  its  second 
argument,  then 

bC 

- =/2(w,w).  (3) 

Similarly,  the  derivative  of  phase  velocity  with  respect  to  to  is  given  by 

^ =/,(co,  m).  (4) 

bU\ 

Now  let  U and  be  values  of  the  functions#  and  j?, 

U = #(o>,  m)  (5) 

bill 

= g2(w,m).  (6) 


I 


and  C which  in 
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is  derived  from  the  relationship  between  U 

'Cl 


c't/l  , iC 

The  relationship  between  and 


i m 


c 

ui  i c 

C i (ui 


<HC. 


. in)- 


(7) 


From  (7),  gUu.  m)  is  recognized  to  be 

y(oj,  m)  - </>( /(io.  m),/,(u>,  m),  <o) . (8) 

Application  of  the  chain  rule  results  in 

8:  = 0i/j  + 0j/ij  (9) 

where  the  function  /, , is  the  derivative  of/,  with  respect  to  its  second  argument.  Since 
/,  and  /',  are  differentiable. 


C /r'Cl  \|  „ (CC\  \J 

; a/ii(«>,  m)  = /2,(w.  m)  s ( . 

< zn  < «)  \< 

?K) 


Explicit  differentiation  of  </>  yields 

4> , 


f/1 


ID. 


- rl 

Therefore,  using  (101.  (II),  and  (12),  equation  (9)  can  be  rewritten  as 


The  expression 


i'U 

V 

, v\ 

CC 

U2  c 

(CC  \ 

2- 

+ (!>  . 

1 

r‘/H 

„ " f'1 

cm  |,u  C cio 

y/m 

III  i'll 

V cm 


<c  Au /."*1 ) 

r‘m10  C i'w\C  Cm\J  „ 


(10) 

(ID 

(12) 

(I-)) 

(14) 


relating  logarithmic  partial  derivatives  can  be  derived  in  a similar  way. 

Novotny  ( 1970)  derived  an  expression  equivalent  to  ( 1 .7).  He  evaluated/, . analytically 
in  terms  of  derivatives  of  Fin  (2),  but  for  Love  waves  only. 

The  next  section  describes  a technique  for  accurately  approximating  (13). 

Numerical  Approximations 

This  section  describes  an  algorithm  which,  for  a given  frequency,  requires  the  solution 
for  only  two  roots  of  equation  (2)  to  get  the  group- velocity  partials  w ith  respect  to  all  of 

the  model  parameters.  The  method  described  is  to  approximate  'jj 


CC 


numerically  differentiating 


Cm 


Standard  Thomson-Haskell  matrix  calculations  are 


used  to  gel  C,  V,  and 


cm 


. Since  equations  ( 1 3)  and  (14)  depend  only  on  the  relationship 


between  U and  C in  (7),  clearly  this  method  can  be  applied  to  any  Rayleigh-  or  Love-wave 
mode. 
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Suppose  the  quantities 


. and 


U 


»,i  = and  *’>  i - m0e  i.  If  the  quantities  (', , U, , and 


are  desired  at  the  frequency  w0.  Define 
C 


for  each  model 


parameter  are  computed  at  m, . and  if  C , , L , . and 


. C 


are  computed  at  w , . then 


the  follow  mj!  approximations  can  be  assigned  to  the  central  frequency  <u„ 

C0  = I 2 (C.+C  ,) 

C„  = I 2 (C.  + IC,) 

<CJ 


■ (n 

t m i 


i-C 

= 12  1 

+ 1 

V"” 

cm 

_ <«(s  1 f'o 

.0  <"o  V C„/  < m 


v°'| 

m v 0 | 

'<ic, 

_<c_, 

\ 

an 

.0  <"l 

- 

\ 

2d 

/ 

(15a) 

(15b) 

(15c) 

(I5d) 


A central  difference  is  used  in  (15d)  because  this  is  the  most  accurate  two-point 


approximation  tow 


<i»\,m 


Using  the  averages  ( 1 5a)  to  ( 1 5c)  avoids  the  calculation 


of  a root  at  <o0.  Also,  the  first  of  the  two  roots  calculated  provides  an  excellent  initial 
guess  for  the  second  root,  which  speeds  up  the  calculation;  given  C , and  V.  ,.  a very 
accurate  guess  to  C,  is  C.  ,[l  + 2<S(!  — (C_ , 'U  ,))). 

The  counterparts  to  equations  (15c)  and  (I5d)  for  the  logarithmic  partials  in  (14)  arc 
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Choosing  d.  The  error  in  the  approximation  (I5d)  can  be  decomposed  into  two  parts; 
e = c,  e,  arises  from  the  nonlinearity  of  C,  U,  and 


i in 


as  a function  of  log  w and 


can  be  shown  to  be  on  the  order  of  d1;  e2  is  the  result  of  round-off  errors  in  C , , Ut, 


< c.l 

cC  , 

# . <‘CY 

, 'C  , 

.c  .. 

C , and  ' 

and  behaves  as  d .It  the  errors  in 

and 
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(17) 


r,  and  i , . respectively,  then  for  small  d 

Uo*  ( «i  -®-i\ 

‘ J ' C„J  [ 2d  ) ■ 

A knowledge  of  r,  and  r.  , enables  a lower  limit  on  d to  be  established  to  ensure  that 
t>2  be  within  a specifievl  tolerance.  If  an  upper  limit  on  |c,|  and  |e.,|,  say  can  be 
determined  based  on  the  precision  of  the  phase-velocity  partial  calculations  performed, 
then  an  upper  limit  on  |e2|  is 


t’j  max  = 


U02  tm«« 

('o'  <5 


(IS) 
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r 

i 


Taking  e2  max  to  be  the  largest  tolerable  |e  ,|,  a criterion  for  choosing  6 is 


rf  > 


e:  max 


(19) 


since  U0  ^ C0.  For  example,  if  * ^ 1 

<m 


if  M decimal  figures  are  desired  for 


and 

r't'o 

'em 


are  accurate  to  N decimal  figures,  and 


I" 


, then  £ma&  = 1/2  10  N,  e2  max  = 1/2  10~M, 


and  (19)  becomes 


5 > I0'*~v.  (20) 

If  M is  sufficiently  less  than  N,  then  <5  will  be  small  enough  to  guarantee  both  a small  e, 
and  a small  e , . The  criterion  in  (20)  can  be  applied  to  all  modes  since  <5  is  a perturbation 
of  log  <0. 


TABLE  1 

CANSD  Canadian  Sun  id  Modi-l 


1 

h,*  (km) 

</,t  (km) 

0,  (km.  sec) 

p,  (gm  cm1) 

1 

6.0 

3.0 

3.47 

2.70 

2 

10.5 

11.25 

3.64 

2.80 

3 

18.7 

25.85 

3.85 

2.85 

4 

80.0 

75.2 

4.72 

3.30 

5 

100.0 

165.2 

4.54 

3.44 

6 

100.0 

265.2 

4.51 

3.53 

7 

80.0 

355.2 

4.76 

3.60 

8 

00 

— 

5.12 

3.76 

* h,  is  layer  thickness, 
t d,  is  depth  to  center  of  ith  layer. 


Examples 


The  algorithm  in  (15),  hereafter  referred  to  as  the  A/*  (delta  partial)  method,  was  tested 
for  fundamental-mode  Love  waves  on  a continental  shield  model  and  for  fundamental- 

( 'C' 

mode  Rayleigh  waves  on  an  ocean  model.  In  both  cases  V,  C,  and  at  the  perturbed 

im  „ 

frequencies(to, andw_,)wcre  calculated  withcomputer  programs  written  by  D.  Harkrider 
for  implementing  methods  described  in  Anderson  (1964),  Takeuchi  et  at.  (1962), 
and  Takeuchi  et  at.  (1964). 

Example  I:  Love  waves.  Using  the  A P method  with  6 = 0.001,  Love-wave  group- 
velocity  partial  derivatives  with  respect  to  shear  velocity  (/))  and  density  (/>)  were  com- 
puted at  periods  of  20  and  40  sec  for  a version  of  the  CANSD  Canadian  shield  model 
(Brune  and  Dorman,  1963),  shown  in  Table  I.  The  same  derivatives  were  computed 
analytically  by  Novotny  (1970).  The  AP  results  are  tabulated  in  Table  2 and  agree  in 
each  case  with  the  five  decimal  figures  given  by  Novotny. 

Example  2:  Rayleigh  waves.  Again  using  equations  (15)  with  <5  = 0.001,  Rayleigh- 
wave  group-velocity  partials  for  a version  of  the  Anderson  ocean  model  (Harkrider. 
1970),  given  in  Table  3,  were  computed  for  seven  periods.  (For  the  purposes  of  this 
example,  the  model  was  taken  as  a flat,  rather  than  spherical,  earth.)  The  group-velocity 
and  phase-velocity  partial  derivatives  with  respect  to  p for  each  layer  at  66,  100,  and 
200  sec  are  shown  in  Figure  I . 
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TABl  I 2 

I’hasi- anii CiKorp-Vi UK  m I’akiiai  Dirivaiimsior  (ANSI)  Mddii 


i 

iC 

<’■*1 

?c 

?pt 

?t/ 

?U 

Period  — 20  sec 1 

C =»  4.00710  kmisec'.  V - 3.5273 

2 km  \ee 

1 

0.25048 

-0  05142 

0.34610 

0.03708 

•> 

0.40458 

0.03721 

0 56217 

0 01013 

3 

0.40666 

0 03642 

0.24470 

0 06446 

4 

0.15636 

0.04255 

-0.17480 

0.03776 

s 

00002: 

0.00005 

0.00205 

0(XX14(. 

6 

0.0 

0.0 

0.0 

0.0 

7 

0.0 

0.0 

0 0 

0.0 

Period  40  sec'. 

C - 4.40204  kmisec;  ( 4.01515  Am  tiv 

1 

0.04234 

-0.03562 

0.26627 

-0.08645 

■> 

0.16802 

-0.04354 

0.47105 

- 004121 

3 

0.26171 

- 0.02070 

0.64765 

0.01354 

4 

0.48448 

0.07374 

0.37427 

0.14447 

5 

0.10746 

0.0081.3 

-0.28054 

— 0.00354 

6 

0.01742 

0.00136 

0.10107 

-0.00573 

7 

0.00147 

0.01X127 

-0.01280 

0,00223 

TABl  1 3 

Anui  rsun  (Vt  an  Mom  i 

i A, 

4, 

8, 

Pi 

0 5. (XXXI 

2.50 

1 5200 

0.0 

1 03(H) 

1 1 (XXX) 

5.50 

2.1000 

1 (XXX) 

2. 1000 

2 5.0000 

8.50 

6.4100 

3.7000 

3.0660 

3 4,  (XXX) 

15.50 

8.1100 

4 6060 

3.4000 

4 5. (XXX) 

22.50 

8.1200 

4.6110 

3.4000 

5 1 5. (XXX) 

32.50 

8.1200 

4.61(H) 

3.4000 

6 10  0000 

45.00 

8.0100 

4.5600 

3.3700 

7 100000 

55.00 

8.0100 

4.5600 

3.37(H) 

8 10.0000 

65.00 

7.4500 

4.5600 

3.3700 

4 lO.(XXX) 

75.  (X) 

7.4500 

4.56(H) 

3.3700 

10  10.0000 

85. (X) 

7.7 100 

4.4000 

3.37(H) 

II  10  (XXX) 

45.00 

7.7 100 

4.4000 

3.37(H) 

12  (O.(KKK) 

105.00 

7.6 H00 

4. 34(X> 

3.3300 

13  10.0000 

1 1 5.00 

7.6800 

4.3400 

3.33(H) 

14  10.0000 

1 25.00 

7.7770 

4.3400 

3.3300 

15  10  0000 

135.00 

7.7770 

4.34(H) 

3.3300 

16  10.0000 

145.00 

7.8500 

4.34(H) 

3.3300 

17  10.0000 

1 55.00 

7.8500 

4.3400 

3.3300 

18  10.0000 

165.00 

8.1  (XX) 

4.45(H) 

3.3300 

14  K).(XXX) 

175.00 

8.1  (XX) 

4.45(H) 

3.3300 

20  10.0000 

185.00 

8.1200 

4.4500 

3.3300 

21  10.0000 

145.00 

8.1200 

4.4500 

3.3300 

22  10.0000 

205.  (X) 

8.1200 

4.45(H) 

3.3300 

23  10.0000 

215.00 

8.12(H) 

4.45(H) 

3.3300 

24  10.0000 

225.00 

8.12(H) 

4.45(H) 

3.3300 

25  10.0000 

235.00 

8.12(H) 

4.45(H) 

3.3300 

26  I0.(XXX> 

245.00 

8.1200 

4.4500 

3.3300 

27  10.0000 

255.00 

8.12(H) 

4.4500 

3.3300 

( onlmued 
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TABLF.  3 — Ctmlinufil 


i 

h. 

</, 

*1 

0, 

<6 

28 

10.0000 

265.00 

8.1200 

4.4500 

3.3500 

29 

10.0000 

275.00 

8 1200 

4 4500 

3.3500 

30 

10.0000 

285.00 

8.1200 

4.4500 

3.3600 

31 

10.0000 

295.00 

8.1200 

4.4500 

3.3600 

32 

10.0000 

305.00 

8.1200 

4.4500 

.3.3700 

33 

10.0000 

315.00 

8.1200 

4 4500 

.3.3700 

34 

10.0000 

325.00 

8.1200 

4.4500 

3.3800 

35 

10. 0000 

335.00 

8.1200 

4 4500 

3.3800 

36 

10  0000 

345.00 

8.2400 

4 5000 

3.3900 

37 

10.0000 

355.00 

8.2400 

4.5000 

3.3900 

38 

10,0000 

365.00 

8.3000 

4.5300 

3 4400 

39 

10.0000 

375.00 

8.3600 

4.5600 

3.5000 

40 

1 0.0000 

385.00 

8.3600 

4.5600 

3.5000 

41 

lO.(XKK) 

395.00 

8.7500 

4.7950 

.3.6840 

42 

1 5.0000 

407  50 

8.7500 

4.7950 

3.6840 

43 

20.0000 

425.00 

9.1500 

5.04(8) 

3.8800 

44 

10.  (XXX) 

440.00 

9.4300 

5.2170 

3.9000 

45 

20.0000 

455.00 

9.7600 

5.4000 

3.9200 

46 

25.0000 

477.50 

9.7650 

5.4000 

3.9330 

47 

25.0000 

502.50 

9.7750 

5,4000 

.3.9480 

48 

25.0000 

527.50 

9.7800 

5.4000 

.3.9600 

49 

25.0000 

552.50 

9.7840 

5.4(XX) 

.3.9880 

50 

25.0000 

577.50 

9.7880 

5.4000 

4.0220 

51 

25.0000 

(8)2.50 

9.7920 

5.4000 

4.0560 

52 

25.0000 

627.50 

9.7960 

5.4000 

4.0900 

53 

25.0000 

652.50 

9.8000 

5 4000 

4.1200 

54 

25.0000 

677.50 

10.16.30 

5.  (XXX) 

4.1650 

55 

25.0000 

702.50 

10.4880 

5.8000 

4.2120 

56 

25.0000 

727.50 

10.8180 

6.1  (XX) 

4.2570 

57 

25.0000 

752.50 

11.1200 

6.2000 

4.3000 

58 

25.0000 

777.50 

11.1350 

6.2050 

4 4750 

59 

25.0000 

802.50 

11.1500 

6.2100 

4.6330 

60 

25.0000 

827.50 

11.1650 

6.2180 

4.7970 

61 

25.0000 

852.50 

11.1809 

6.2.300 

4.9400 

62 

25. (XXX) 

877.50 

11.2240 

6.2500 

4.9425 

63 

25.0000 

902.50 

11.2670 

6.2750 

4.9450 

64 

25.0000 

927.50 

II.  .3 100 

6.2970 

4.9475 

65 

25.0000 

952.50 

11.3500 

6.3220 

4.9500 

66 

25.0000 

977.50 

11.3920 

6.. 3400 

4.9517 

67 

25.0000 

1002.50 

11.4.340 

6.3600 

4.9534 

68 

25.0000 

1027.50 

1 1 .4760 

6.3750 

4.9550 

69 

25.0000 

1052.50 

11.5180 

6.3900 

4.9567 

70 

25. (XXX) 

1077.50 

11.5600 

6.4050 

4.9584 

71 

00 

— 

11. 6000 

6.4210 

4.9600 

The  partials  with  respect  to  /i  in  six  of  the  model's  layers  were  also  computed  by  the 
•'brute  force"  approximation  for  live  periods.  A central  difference  scheme 

(W 


M 

2A/f 


(2\) 


was  used,  where  A V is  the  difference  between  the  group  velocities  for  models  with  /f 
perturbed  by  + A/I  and  - A/f  in  a given  layer.  A A /i  of  0.1  was  used  for  each  ol  the  6 
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layers.  The  accuracy  of  the  l partials  computed  this  way  is  difficult  to  judge,  but  repeating 
the  layer  6 calculation  with  A/#  = 0.05  revealed  that  the  errors  were  on  the  order 
of 0.00002. 

It  can  be  seen  from  Table  4 that  the  agreement  between  the  brute  force  and  AP  group- 
velocity  partials  is  very  close.  For  layer  6 the  agreement  is  closer  when  A/i  = 0.05  than 
when  A p = 0.1.  Therefore,  we  infer  errors  of  0.00002  or  less  in  the  AP  answers. 


Discussion 

The  AP  method  is  an  efficient  way  to  compute  group-velocity  partial  derivatives  with 
an  accuracy  sufficient  for  use  in  inversion  calculations.  This  capability  increases  the 
feasibility  of  using  available  group-velocity  data,  either  separately  or  jointly  with  phase- 
velocity  data,  in  the  determination  of  earth  structure  by  linear  inverse  methods. 

PARTIAL  DERIVATIVES  (ito-3) 


Fio.  I.  Phase-  and  group- velocity  partial  derivatives  with  respect  to  fi  at  periods  of  hh.  100,  and 
200 sec  for  Anderson  ocean  model.  Partial  derivatives  are  normalized  bv  layer  thickness;  i.c.,  (I  A,) 
(iC/ffii)  and  t.l/A,K3t//c,(J()- 

It  is  very  desirable  to  be  able  to  invert  group-velocity  data  because  of  the  extensive 
areal  coverage  of  existing  group-velocity  observations  as  compared  to  phase-velocity 
observations,  particularly  at  periods  which  sample  the  lithosphere.  Also,  a very  large 
percentage  of  available  higher-mode  dispersion  data  consists  of  group-velocity  measure- 
ments. Moreover,  group-velocity  determinations  commonly  can  be  more  accurate 
than  phase-velocity  determinations.  This  is  particularly  true  if  a single  source-station 
pair  is  used  or  if  Fourier  spectral  methods  involving  long  time  windows  are  required  for 
phase-velocity  estimates  while  narrow-band  filters  suffice  for  group  velocity. 

In  addition  to  their  use  in  finding  models  consistent  with  group-velocity  observations, 
partial  derivatives  can  be  used  to  generate  averaging  kernels.  With  averaging  kernels, 
experiments  can  be  designed  to  determine  the  most  favorable  frequencies  to  observe  and 
the  accuracy  of  group-velocity  data  needed  to  resolve  radial  or  lateral  variations  in  earth 
structure. 


CALCULATION  OF  RAYLEIGH-  AND  LOVE-MODE  GROl'P-VELOCITY  PAR  M AI 


TABLE  4 

Comparison  of  AP  ani>  Brutf  For.  t (B  f .)  Rayiik.m- 
wavf  Group-Velocity  Pariials  with  Rispict  to  ft 
|(?l//'70U  in  Six  I ayers  o(  Anderson  Ocf.an  Model* 


I B F. 


Period  = I40.vtv;  ('  = 4.1782  km,  see, 

U — 3.7468  km  sec 

6 

0.01378  (A/I  = 0 051 

A 

0.01379 

0 01378 

II 

0 05285 

0.03281 

IA 

O.OaOaa 

0.06063 

21 

0.04024 

0 0402 ' 

31 

-0.00058 

0.00058 

41 

-0.01140 

-0.01139 

Period  --  100  sec:  C = 4 0751  A 
l = 3.9087  km  see 

m see : 

A 

0.03559  (A ft  = 0.05) 

A 

0.035A0 

0.03557 

II 

0.08448 

0.08444 

IA 

0.05468 

0.03466 

21 

0.01386 

0.01387 

31 

-0  01480 

0.01480 

41 

-0.00770 

0.00769 

Period  *=  10  sec;  C = 4.0406  km  see. 

V = 4.0003  km  sec 

A 

0.08A  1 7 (A/I  = 0.0ft 

A 

0.08A22 

0.0861 5 

II 

0 09409 

0.09407 

IA 

0.01697 

0 01 698 

21 

-0.01484 

0.01483 

31 

-0.01113 

0.01113 

41 

-0.00207 

-0.00207 

Period  = 33  sec;  C - 4.0.V50  km,  sec; 

U = 4.0051  km  sec 

A 

0.23231  (Aft  = 0.05) 

A 

0.23234 

0.23231 

II 

- 0.02903 

-0  02903 

IA 

- 0.02447 

0.02446 

21 

— 0.00642 

0.00642 

31 

-0.00023 

— 0.00023 

41 

— O.tXXXX) 

-0.00000 

Perio,/  = 1 A. 5 see-.  C 3.8375  k 
V = 2.9741  Am  sec 

m see  . 

A 

0.1087a  (Aft  0.05) 

A 

-0.10877 

-0  10876 

II 

— 0.00783 

— 0.00782 

IA 

-0.00013 

-0.00013 

21 

O.tXXXX) 

0.00000 

31 

0.00000 

0 tXXXX) 

41 

0.00000 

O.tXXXX) 

• The  brute  force  results  were  computed  with  AN  0 I 
except  where  0.05  is  indicated. 
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An  inversion  program  is  currently  being  developed  that  incorporates  the  AP  method 
of  computing  partial  derivatives  and  allows  inversion  of  group  and  phase  velocity 
separately  or  in  combination  for  both  Rayleigh  and  Love  modes. 
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